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Abstract 

A  technique  is  developed  to  determine  the  orbit  of  a  sunlight  illuminated  satel¬ 
lite  passing  through  the  field-of-view  of  a  sensor  platform  in  a  Highly  Elliptical  Orbit 
(HEO)  and  Geosynchronous  orbit  (GEO).  The  technique  develops  two  different  meth¬ 
ods  of  initial  orbit  determination.  The  first  relies  on  the  Gauss  initial  orbit  determi¬ 
nation  method  to  develop  an  estimate  of  the  state  from  angular  data.  The  second 
method  relies  on  positional  data  of  the  target  relative  to  the  Earth’s  background  to 
determine  an  estimate  of  the  state.  These  estimates  are  then  refined  in  a  non-linear 
least  squares  routine.  This  estimate  of  the  state  is  then  used  to  identify  the  target 
from  the  Air  Force  Space  Command  satellite  catalog. 

It  was  found  that  the  Gaussian  initial  orbit  determination  method  produced 
reasonable  estimates  of  the  state  given  data  sets  larger  than  30  seconds.  The  second 
method  proved  to  be  useful  when  only  very  limited  data  sets  were  available  and  had 
no  limitations  on  data  size.  The  least  squares  process  was  able  to  compensate  for 
errors  of  as  much  as  5  degrees  in  the  initial  estimate  of  the  state  values.  In  all  HEO 
test  cases  the  estimate  of  the  state  vector  provided  enough  information  to  identify  the 
target  satellite  from  the  satellite  catalog.  In  GEO  test  cases  it  proved  impossible  to 
identify  other  objects  in  GEO  given  the  limited  data  available.  The  method  developed 
also  proved  unreliable  in  identifying  orbits  with  eccentricities  greater  than  .35.  Over 
82  percent  of  the  objects  listed  in  the  Air  Force  Space  Command  satellite  catalog  can 
be  successfully  identified  in  orbit  using  this  method. 
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Space  Based  Satellite  Tracking  and  Characterization 


Utilizing  Non-Imaging  Passive  Sensors 

I.  Introduction 

1.1  Objective  and  Scope 

Satellite  based  sensors  looking  down  at  the  Earth’s  surface  occasionally  observe 
reflected  light  from  an  object  passing  through  the  image  which  is  moving  too  fast 
relative  to  the  background  of  the  image  to  be  located  within  the  atmosphere.  These 
objects  are  commonly  called  fastwalkers.  This  term  refers  to  any  orbital  object  seen 
passing  through  the  held  of  view  of  an  Earth  observing  sensor  which  is  suspected  of 
being  in  orbit.  The  objective  of  this  thesis  is  to  develop  a  method  to  determine  the 
orbit  of  these  objects  from  a  sensor  in  a  highly  elliptical  orbit  (HEO)  and  geosyn¬ 
chronous  orbit  (GEO)  and  then  match  them  to  an  object  in  the  Air  Force  Space 
Command  (AFSPC)  catalog  of  satellites. 

In  order  to  conduct  this  analysis  a  few  basic  assumptions  must  be  made.  First, 
the  sensor  ephemerids  are  known  exactly  since  the  position  of  the  observer  is  necessary 
to  relate  the  observations  to  a  location  in  space.  Second,  the  targets  are  assumed  to 
be  non-thrusting  bodies  in  low  Earth  orbit  (LEO),  meaning  they  are  at  or  near  the 
Earth’s  surface  below  2000  km  in  altitude.  This  assumption  is  reasonable  for  a  HEO 
satellite  because  the  number  of  objects  in  a  HEO  orbit  is  very  small,  so  very  few  non- 
LEO  spacecraft  are  likely  to  appear  in  the  HEO  sensors  held  of  view.  This  assumption 
is  not  valid  for  a  GEO  sensor.  Therefore,  we  will  need  to  differentiate  between  targets 
in  LEO  and  GEO  when  using  observations  from  a  GEO  platform.  Lastly,  since  the 
objects  are  in  a  LEO  orbit  and  the  sensor  held  of  view  is  limited,  we  will  assume 
that  the  maximum  amount  of  time  for  which  a  data  track  is  visible  is  less  than  4 
minutes.  Given  this  brief  period  of  observational  data,  it  will  be  assumed  that  orbital 
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perturbations  will  have  no  visible  effect  on  the  orbit.  For  this  reason  the  two  body 
equations  of  motion  are  sufficient  to  estimate  the  orbital  element  set. 

The  simulated  data  for  this  project  was  generated  by  the  Satellite  Tool  Kit 
(STK)  available  from  Analytical  Graphics  Incorporated  (AGI).  This  software  is  cur¬ 
rently  among  the  most  versatile  available  and  is  commonly  used  for  orbital  modeling 
within  the  Army  and  Air  Force  Space  Community.  Version  7  of  STK  was  used  in  the 
generation  of  all  orbital  data. 

1.2  Research 

In  his  thesis  Osedacz  [9]  states  that  a  database  of  fastwalker  observations  exists 
with  data  from  as  far  back  as  1972.  He  uses  this  data  from  obtained  ballistic  missile 
early  warning  satellites  in  GEO  to  attempt  to  determine  the  orbit  of  these  fastwalkers. 
His  conclusion  was  that  determining  the  orbit  of  other  fastwalkers  in  GEO  was  im¬ 
possible.  This  supported  the  conclusions  of  an  earlier  work  conducted  by  Lem  Wong 
on  the  inability  to  determine  the  orbit  of  these  objects.  The  research  conducted  by 
Wong  is  no  longer  available,  but  according  to  Osedacz,  the  problem  was  found  to  be 
unsolvable  in  a  GEO  orbit.  This  conclusion  stemmed  from  the  fact  that  most  of  the 
fastwalkers  were  determined  to  also  be  in  Geosynchronous  orbit  with  the  sensor.  The 
result  was  an  inability  to  determine  the  range  to  the  target. 

The  two  previously  mentioned  works  dealing  with  fastwalkers  focused  on  sensors 
in  GEO  orbit.  In  a  recent  paper  published  by  Li,  Guo,  and  Zhou  at  the  Chinese 
National  University  of  Defense  Technology,  a  Kalman  filter  was  used  to  successfully 
determine  an  estimate  of  the  orbital  elements  of  an  object  visible  from  a  sensor  in 
HEO  orbit  [7].  The  Kalman  filter  needed  over  30  minutes  of  tracking  data  to  converge 
on  a  solution  that  was  accurate  to  within  4  kilometers.  This  method  was  inapplicable 
to  the  central  problem  of  this  thesis,  since  a  sensor  in  HEO  orbit  cannot  normally 
observe  an  object  in  LEO  orbit  for  that  amount  of  time.  Very  little  research  exists 
using  angles-only  data  from  a  space  platform  and  almost  nothing  on  the  problem  of 
fastwalker  orbit  determination. 


2 


1.3  Approach 

In  addition  to  the  sensor  ephemerids,  it  was  assumed  that  the  time  of  each 
observation,  the  azimuth  and  elevation  to  the  target  in  the  sensor  frame,  and  the 
latitude  and  longitude  of  the  target  against  the  background  of  the  Earth  were  available 
data.  The  general  method  of  correlating  observed  tracks  to  entries  in  the  Air  Force 
Space  Command  catalog  was  to  use  one  of  two  initial  orbit  determination  techniques 
and  then  to  use  a  least  squares  process  to  refine  the  estimate  of  the  state  vector.  The 
catalog  could  then  be  screened  for  possible  matches  by  successively  applying  filters 
based  on  the  estimates  of  the  elements  in  the  state  vector. 

The  state  vector  consisted  of  four  elements.  Three  of  which  are  necessary  to 
identify  the  target.  The  elements  of  the  state  vector  are:  right  ascension  of  the 
ascending  node  (RAAN),  inclination  of  the  orbit  plane,  an  argument  of  latitude  which 
is  equal  to  the  mean  anomaly  of  a  circular  orbit  and  an  estimate  of  the  altitude  of  the 
target  fastwalker.  For  the  purposes  of  this  research,  both  the  node  and  the  inclination 
were  treated  as  constants.  They  were  then  used  to  search  the  catalog,  assuming  that 
the  most  recent  catalog  is  used  so  that  the  change  in  RAAN  is  small. 

The  first  initial  orbit  determination  technique  used  the  latitude  and  the  longi¬ 
tude  of  the  background  image.  Spherical  trigonometry  was  then  used  to  determine 
the  state.  The  second  initial  orbit  determination  technique  was  the  classic  initial  orbit 
determination  technique  developed  by  Gauss.  This  technique  used  the  azimuth  and 
elevation  data  to  estimate  the  state  vector. 

Once  an  initial  estimate  of  the  state  vector  was  obtained,  a  non-linear  least 
squares  process  was  applied  to  refine  the  estimate  using  the  rest  of  the  data  points 
available  in  a  given  track.  This  refined  estimate  of  the  state  was  then  used  to  search 
the  catalog  for  the  right  ascension  of  the  ascending  node  and  inclination  near  the 
estimate  of  the  state.  Once  a  candidate  was  found,  it  could  be  propagated  to  the 
epoch  time  of  the  estimated  state  and  it’s  argument  of  latitude  compared  to  that  of 
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the  estimated  state.  A  match  would  then  allow  us  to  confidently  state  that  we  had 
identified  the  orbital  object. 
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II.  Orbit  Generation  Equations  and  Methods 


This  chapter  presents  the  equations  necessary  to  develop  an  accurate  model 
of  an  orbiting  object  for  the  purpose  of  determining  its  orbit.  The  element  sets, 
coordinate  frames  and  transformations  used  will  be  discussed  in  detail.  This  chapter 
will  also  focus  on  describing  the  problem  geometry  and  develop  two  different  initial 
orbit  determination  techniques.  The  first  initial  orbit  determination  technique  that 
will  be  discussed  is  the  classic  Gauss  method,  which  uses  three  measurements  of 
azimuth  and  elevation  in  the  sensor  frame.  The  second  technique  discussed  uses 
latitude  and  longitude  data  from  the  background  of  the  image  produced  by  the  sensor. 
Descriptions  of  both  techniques  and  their  development  are  discussed  in  this  chapter. 

2.1  Problem  Geometry 

Understanding  the  geometry  of  the  problem  is  fundamental  to  solving  it.  Both 
initial  orbit  determination  methods  rely  on  the  same  fundamental  equation.  If  we 
let  the  vector  r  represent  the  position  of  the  object  with  respect  to  the  center  of  the 
Earth,  p  is  the  range  of  the  target  satellite  from  the  sensor,  and  R  is  the  position  of 
the  sensor  then  we  can  write  by  observation  that 


r  =  p  +  R  =  pL  +  R 


(2.1) 


where  L  is  a  vector  of  direction  cosines,  as  shown  in  Figure  2.1. 

2.2  Element  Sets 

2.2. 1  Classical  Orbital  Elements.  To  describe  an  orbit  under  the  assumption 
of  a  two-body  gravity  model,  six  or  more  parameters  are  required  in  addition  to  an 
epoch  time.  For  example,  given  some  time,  to,  and  the  position  and  velocity  vectors 
of  a  satellite,  we  can  fully  describe  the  current  state  of  a  satellite.  However,  due  to  the 
limited  nature  of  the  data  obtained  from  our  observations  and  the  need  to  limit  the 
number  of  values  in  our  state  vector  when  using  non-linear  least  squares,  the  Classical 
Orbital  Element  set  was  chosen  as  the  parameter  set  in  this  work.  In  this  element  set 
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Sensor 


Figure  2.1:  Problem  Geometry 


five  of  the  six  elements  are  constant  and  one  varies  periodically  through  a  single  orbit 
of  the  observed  spacecraft.  For  this  reason,  the  Classical  Orbital  Elements  simplifies 
the  process  of  developing  an  estimate  of  the  satellite  state.  Although  many  sources 
define  this  orbital  set,  Wiesel  [16]  provides  a  solid  development,  which  we  will  vary 
from  only  slightly  here. 

The  classical  orbital  elements  are: 


a  -  The  Semimajor  Axis 
e  -  Eccentricity 
i  -  Inclination 

fl  -  Right  Ascension  of  the  Ascending  Node 
w  -  Argument  of  Perigee 
v  -  True  Anomaly 

When  describing  the  Classical  Orbital  Elements,  it  is  necessary  to  first  determine  the 
shape  of  orbits  around  a  body.  An  orbit  around  a  body  can  be  described  as  a  conic 
section,  which  is  generated  by  slicing  a  circular  cone  with  a  plane.  Five  types  of  conic 
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Figure  2.2:  The  Classical  Orbital  Elements 


sections  can  be  formed  in  this  way,  but  we  are  only  concerned  with  two,  the  circle  and 
the  ellipse.  The  semi  major  axis,  a,  is  a  constant  which  describes  half  of  the  major 
axis  of  an  ellipse.  The  shape  of  the  ellipse  is  defined  by  the  eccentricity,  e,  which  is  a 
value  in  the  interval  0  <  e  <  1  and  for  a  circle  e  =  0. 

The  inclination  of  an  orbit  refers  to  the  tilt  of  the  orbital  plane.  This  angle  is 
measured  from  the  unit  vector  described  by  the  Z  axis  to  the  angular  momentum 
vector,  h,  and  varies  from  0  <  i  <  180  degrees,  as  shown  in  Figure  2.2.  It  can  be 
defined  in  terms  of  these  quantities  by, 


cos  ( i ) 


Z-h 


(2.2) 


The  right  ascension  of  the  ascending  node,  fl,  is  the  angle  in  the  equatorial 
plane  measured  positively  from  the  X  axis  to  the  point  on  the  orbit  where  the  satellite 
crosses  the  equator  moving  from  South  to  North.  Values  for  the  ascending  node  range 
from  0  <  <  360  degrees.  All  inclined  orbits  also  have  a  descending  node  where  the 
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satellite  crosses  the  equator  moving  from  North  to  South.  The  line  connecting  these 
two  points  is  called  the  Line  of  Nodes. 

The  final  two  elements  of  the  set  are  interrelated.  The  Argument  of  Perigee,  w, 
is  the  angle  formed  between  the  Line  of  Nodes  and  the  point  at  which  the  orbit  is 
closest  to  the  Earth,  called  perigee.  Values  for  the  Argument  of  Perigee  range  from  0 
to  360  degrees.  The  True  Anomaly,  v,  is  undefined  for  perfectly  circular  orbits,  which 
the  least  squares  model  in  this  thesis  assumes.  For  circular  orbits  this  value  is  replaced 
with  a  measurement  of  the  angle  between  the  right  ascension  of  the  ascending  node 
and  the  actual  position  of  the  satellite  at  some  epoch  time.  For  this  thesis,  we  will 
refer  to  this  as  the  Argument  of  Latitude.  For  non-circular  orbits,  the  True  Anomaly 
is  not  measured  from  the  line  of  nodes  but  from  the  perigee  point  as  pictured  in  Figure 
2.2.  Values  for  both  also  range  from  0  to  360  degrees. 


ion  1st  derivative 


Bstar  Epoch  Elem  # 


Card  #  International  Designator  Epoch  YYDDD  Mean  Motion  2nd  derivative 
/  Common  Name  /  Mean  I 

/  l 

l£S  (ZARYA) 

1  25544U  98067A  07270.64846184  .00007688  00000-0  54207-4  0  607 

2  25544  051 .6348  295.2535  0002708  045.3662  054.4413  15.75238918506947 
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Figure  2.3:  Two  Line  Element  Set 

2.2.2  Two  Line  Element  Set.  While  the  Classical  Orbital  Elements  are  used 
widely,  the  Air  Force  Space  Command  catalog  of  satellites  is  published  using  a  set  of 
elements  commonly  referred  to  as  the  Two  Line  Element  set  (TLE).  An  example  of 
a  TLE  is  shown  in  Figure  2.3.  This  system  of  describing  orbits  is  a  relic  of  the  need, 
before  the  modern  age  of  computing,  to  use  punch  cards  in  order  to  have  the  data  read 
by  computers.  Space  Command’s  catalog  is  the  most  up  to  date  and  reliable  catalog 
of  objects  in  orbit  around  the  Earth.  For  this  reason,  it  is  used  as  the  database  for 
our  search. 


2. 3  Coordinate  Frames 


Several  Coordinate  Frames  are  used  for  the  purposes  of  generating  data  and 
describing  the  orbits  of  the  objects  within  them.  The  transformations  between  these 
coordinate  frames  and  their  use  is  described. 

2.3.1  Earth  Centered  Inertial  Frame.  The  Earth  Centered  Inertial  Frame 
(ECI)  originates  at  the  center  of  the  Earth  and  is  designated  in  Figure  2.2  by  the 
letters  XYZ.  The  fundamental  plane  of  this  coordinate  frame  is  the  Earth’s  Equator. 
The  X  axis  points  towards  the  vernal  equinox.  The  Y  axis  is  90  degrees  to  the  East 
of  the  X  axis  in  the  equatorial  plane.  The  Z  axis  extends  through  the  North  Pole. 
Another  set  of  labels  that  are  commonly  used  for  these  are  UK. 

2.3.2  Sensor  Frame.  The  sensor  frame  is  a  slightly  modified  Topocentric 
Horizon  Coordinate  System  (SEZ).  This  system  is  commonly  used  for  ground  based 
telescopes  where  the  local  horizon  forms  the  fundamental  plane.  The  S  axis  points 
South  from  the  observer.  The  E  axis  points  East  and  is  undefined  at  the  poles.  The  Z 
axis  points  outward  from  the  center  of  the  Earth.  In  the  sensor  frame  the  coordinate 
system  is  translated  along  the  radial  direction  from  the  surface  of  the  Earth  to  the 
center  of  the  observing  satellite.  See  Figure  2.1  for  an  example  of  the  coordinate 
frame  with  reference  to  the  problem  geometry.  Measurements  taken  in  this  frame  are 
azimuth  (Az)  and  elevation  (El).  Azimuth  is  measured  positively  clockwise  from  the 
-S  axis  and  elevation  is  measured  negatively  down  from  the  plane  formed  by  the  E 
and  Z  axis  where  the  negative  Z  axis  would  be  -90  degrees.  The  line-of-sight  unit 
vector  to  the  target  in  the  sensor  frame,  ,  can  then  be  defined  as: 


- 1 

t-i 

_ 1 

—cos  ( Az )  cos  (El) 

Lhy 

= 

sin  (Az)  cos  (El) 

1 

sin  (El) 
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Figure  2.4:  Right  Ascension  and  Declination 

2.3.3  Right  Ascension  and  Declination  Frame.  When  using  the  Gauss  Initial 
Orbit  Determination  Method  the  measurements  are  taken  as  azimuth  and  elevation  in 
the  sensor  frame  and  need  to  be  transformed  into  right  ascension  and  declination.  In 
order  to  transform  from  the  sensor  frame  to  the  right  ascension-declination  frame,  the 
sensor  frame  unit  vector  found  from  equation  (2.3)  is  needed.  Before  transformation 
between  the  two  coordinate  frames  is  possible,  the  two  angles  described  by  A  and  0 
in  Figure  2.4  need  to  be  found.  Recall  that  we  assumed  we  know  the  position  of  the 
sensor,  R.  R’s  ECI  components  can  be  defined  as  R,t ,  Rj,  Rk ■  Then  a  vector  Rij  may 
be  defined  as 


R, 

Rj 

0 


(2.4) 
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The  angle  0  can  then  be  found  by: 


(2.5) 


We  can  then  similarly  use  either  the  /  axis  unit  vector  or  the  i  component  of  R  to 
find  A  by: 


A  =  acos 


(  Rij 

X 

\ 

Rij 

X 

) 

(2.6) 


Care  must  be  taken  when  using  this  equation,  it  requires  a  quadrant  check  of  the 
position  of  the  sensor  to  determine  if  6  is  greater  than  180  degrees.  With  6  now 
defined  it  is  then  possible  to  transform  the  sensor  frame  measurement  unit  vector 
into  the  right  ascension-declination  frame. 


Li 

sin  (0)  cos  (A)  —sin  (A)  cos  (0)  cos  (A) 

Lhi 

Lj 

= 

sin  (0)  sin  (A)  cos  (A)  cos  (0)  sin  (A) 

Lhj 

Lk 

—cos  (0)  0  sin  (0) 

Lhk 

2.3.4  Geocentric  Latitude  and  Longitude.  Measurements  of  latitude  and 
longitude  can  be  found  using  the  background  of  images  taken  by  the  sensor  platform. 
These  measurements  can  then  be  used  to  approximate  the  position  of  a  low  earth  orbit 
satellite’s  position  near  the  Earth’s  surface.  These  measurements  are  typically  given  in 
terms  of  surface  coordinates  on  the  Earth’s  surface,  referred  to  as  geodetic  coordinates. 
Since  the  Earth  is  not  perfectly  spherical  and  does  in  fact  have  a  slightly  elliptical 
shape  it  is  necessary  to  convert  these  measurements  to  geocentric  coordinates.  First 
we  need  to  define  a  flattening  constant,  Fearth  which  describes  the  elliptical  shape. 


Fearth 


298.257223563 


(2.8) 
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It  is  then  possible  to  find  the  geocentric  latitude  given  the  geodetic  latitude. 

Lat geocentric  QitCLTl  ^(1  dearth  )  tUTl  (-^Otgeodedjc))  (2-9) 

2.4  Initial  Orbit  Determination  Techniques 

The  nature  of  the  data  provides  no  direct  measurement  of  range.  We  are  then 
limited  to  angles  only  initial  orbit  determination  techniques  by  this.  These  techniques 
have  classic  roots  with  the  methods  of  Gauss  and  Laplace  with  more  recent  work  in 
the  held  having  been  done  by  R.H.  Gooding  [4,5]  and  in  a  published  technical  paper 
by  the  National  Aeronautics  and  Space  Administration  (NASA)  [2],  Despite  this,  the 
number  of  unique  methods  developed  since  Gauss  and  Laplace  is  extremely  limited. 
For  over  150  years,  there  were  many  refinements  of  these  two  unique  techniques  devel¬ 
oped.  However,  the  next  significant  work  was  not  published  until  1959  by  R.E.  Briggs 
and  J.W.  Slowley  [1],  although  mention  is  made  by  both  B.G.  Marsden  [8]  and  L.K 
Kristensen  [6]  of  a  little  known  technique  developed  in  1939  by  Vaisala  for  which  I 
could  find  only  broad  descriptions.  The  Briggs  and  Slowley  technique,  relying  on  the 
recent  development  of  computers,  introduced  an  iterative  method  of  angles  only  initial 
orbit  determination.  This  method  formed  the  basis  for  an  entirely  new  method  by  P. 
Escobal  [3]  called  double-r  iteration.  Most  recently,  Vallado  [15]  reviewed  the  three 
primary  techniques  of  Gauss,  Laplace,  and  Escobal  iteration,  recommending  them  as 
the  current  standard  in  initial  orbit  determination.  This  work  considered  only  these 
three  methods  and  R.H.  Gooding’s  recent  work  in  selecting  a  suitable  initial  orbit 
determination  method. 

These  four  techniques  represent  a  broad  range  of  approaches  to  solving  the  initial 
orbit  determination  problem,  but  we  will  focus  only  on  their  effectiveness  in  terms 
of  the  number  and  angular  separation  of  the  observations  that  is  necessary  for  them 
to  function.  The  Gooding  technique  is  an  iterative  technique  developed  specifically 
for  Earth  orbiting  satellites  which  relies  on  angular  data  developed  over  multiple 
revolutions.  For  our  purposes,  this  method  is  impractical  for  the  problem  considered 
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in  this  thesis.  A  thorough  development  of  the  remaining  techniques  is  developed 
by  Escobal  [3].  He  states  that  the  double-r  iteration  technique  is  preferable  but  it 
requires  data  spanning  large  arcs  or  a  number  of  revolutions.  The  standard  Laplacian 
technique  is  a  classic  technique  optimized  for  determining  the  orbits  of  very  distant 
heliocentric  planetary  orbits  [3].  The  Laplacian  method  is  not  considered  robust 
or  reliable  enough  for  satellite  observations  [15].  Escobal  states  that  the  preferred 
method  for  short  arcs  is  Gauss  [3].  However,  this  recommendation  is  not  unanimous. 
A  vigorous  condemnation  was  made  of  the  Gauss  technique  by  L.G.  Taff  [13],  in  which 
he  cited  its  poor  radius  of  convergence  based  on  it’s  use  of  series  functions  but  it  has 
also  been  equally  well  defended  by  B.G.  Marsden  [8]. 

Two  primary  initial  orbit  determination  techniques  will  be  used.  The  first  is  the 
classic  Gauss  Technique  which  relies  on  three  position  vectors  and  their  associated 
times.  The  Gibbs  or  Herrick-Gibbs  technique  is  then  used  to  determine  the  velocity 
vector  at  the  central  position  and  obtain  a  complete  description  of  the  orbit.  The 
Gauss  method  can  be  modified  to  use  either  azimuth  and  elevation  measurements  in 
the  sensor  frame  or  geocentric  latitude  and  longitude  measurements  obtained  from 
the  background  of  the  image  generated  by  the  sensor.  The  second  technique  is  less 
sophisticated  and  uses  the  same  latitude  and  longitude  measurements  applied  on  a 
sphere  to  obtain  a  partial  description  of  the  orbit. 

2-4-1  Gauss  Angles-Only  Initial  Orbit  Determination.  This  development 
and  description  of  the  Gauss  technique  follows  that  of  Escobal  [3]  and  Pressing  and 
Conway  [10].  Since  the  motion  of  the  object  occurs  in  a  plane,  it  is  possible  to  express 
the  position  vector  at  any  time  as  a  linear  combination  of  the  position  vectors  at  any 
two  other  times.  Therefore,  it  should  be  possible  to  determine  a  set  of  scalars  a,  b,  c 
not  all  zero  such  that 


afi  +  fer2  +  cr3  =  0 


(2.10) 
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If  we  solve  for  the  central  date  corresponding  to  r2  and  denote  the  new  coefficients 
by  ci  and  C3, 


r2  =  Cifi  +  c3r3 


(2.11) 


Expressions  for  the  coefficients  in  the  above  equation  can  be  found  by  taking 
the  cross  products  of  the  equation  with  the  first  and  third  position  vectors. 

r2  x  r  3 
r 1  x  r  3 


(2.12) 


c3  = 


\r2  x  rd 


(2.13) 


|r3  x  fi| 

To  solve  for  ci  and  c3,  the  /  series  function  and  g  series  function  are  used.  These 
functions  allow  us  to  express  the  position  vector  at  some  time,  r),  in  terms  of  the 
velocity,  v,  and  position,  ri±j  at  some  other  time. 


r,  =  /i?2  +  Sift  i  r3  =  f3f2  +  g3v2 


(2.14) 


where  /  and  g  are: 


f  _  n  —  t2)2 

h  —  1 - x - h-2  >  93  —  [h  —  h)  ~ 


(£3  —  h  y 
6 


(2.15) 


f  _i  (ti~t2)2u  _  (  +  ^  (h  —  h)3 

/ i  —  l - x - h 2  ,  g  1  —  (t  1  —  t2) - - - 


2  “  ’  J  v  6 
where  we  can  define  h2  using  Earth’s  gravitational  parameter  p, 


(2.16) 


/l2  =  — 3 

|r2|3 


(2.17) 
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Now  we  can  define  a  new  time  variable,  r,  to  simplify  the  time  differences  in  equations 
(2.15)  and  (2.16). 


=  £3  —  t2  ;  T'2  =  t:>,  t\  ;  r3  =  t2  —  t\ 


(2.18) 


If  we  combine  equation  (2.12)  and  equation  (2.14)  we  get; 


Similarly: 


\r2  X  f3 1  \r2  x  (/3r2  +  g3u2)| 

\n  x  r3 1  |(/if2  +  0iu2)  x  (/3r2  +  03v2)| 


Cl  = 


93 

(fi93  —  hgi) 


(2.19) 


t3 

c3  «  — 
r2 


1  + 


lr2  ~  r3  . 
6 


ho 


(2.20) 


Combining  equation  (2.1)  with  equation  (2.11)  we  obtain: 


p2L2  —  C1P1L1  —  c3p3L3  —  cyr sitel  +  c3rSjte3  —  rSjte2 


(2.21) 


where: 

•  ci  and  c3  are  functions  of  |r2|  and  ti,t2,f3. 

•  Li,L2,L3  are  observationally  determined  from  either  the  background  latitude 
and  longitude  or  the  observed  azimuth  and  elevation  angles  from  the  observer 
to  the  target. 

•  rSitei,rsite2-,rSite3  are  the  position  vectors  of  the  sensor  platform  with  respect  to 
the  center  of  the  Earth  and  are  known. 

Equation  (2.21)  represents  a  set  of  three  scalar  equations  with  the  four  unknowns 
P11P21P3  and  |r2|.  Therefore  one  more  independent  equation  must  be  added  to  create 
a  system  that  can  be  solved. 
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(2.22) 


|r2|2  =  pi  +  \rsite2\2  +  2p2-^2  •  rsite2 
We  can  now  rewrite  equation  (2.11)  as: 


Cifi  +  c2f2  +  c3r3  =  0 


(2.23) 


where  c2=-l  or  in  terms  of  equation  (2.21), 


C2P2L2  +  C1P1L1  +  C3p3L3  —  —C\f Slte\  —  C3rsite3  —  C2Tsite2  =  M  (2.24) 


expanded  in  matrix  form, 


L\x  L/2x 

Cl  Pl 

Clpl 

L\ y  L2  y  y 

C2P2 

=  M  =  [L\ 

C2P2 

1 - 

to 

tr-* 

CO 

isj 

C3p3 

C3p3 

(2.25) 


Therefore,  we  can  separate  the  matrix  of  line-of-sight  unit  vectors  by  inverting  it, 


ClPl 

C2P2 

c3p3 


[L\~l  M 


(2.26) 


where 


[L]'1  =  [4  = 


an 

012 

013 

^21 

O22 

023 

031 

032 

033 

(2.27) 


At  this  point  we  need  to  redefine  (2.19)  and  (2.20)  in  terms  of  new  variables, 


Cl 


9  3  _  n 

/1.9.3  -  f3gi  t2 


A\  +  B3h2 


(2.28) 
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T3 

C3  ~  — 
T2 


\  |  frl-Tl)  ' 
6  2 


—  ^4  3  +  -B3^2 


(2.29) 


where  we  then  define  the  vectors  A  and  I?  as: 


iP’  = 


Ai  —1  A3 


;  bt  = 


B\  0  B3 


(2.30) 


We  can  also  define  the  site  vector,  rSite,  at  all  three  times  such  that, 


n  T 


x  ;  y 


^  site  lx 

^  site  2x 

site3x 

T'sitely 

T site2y 

site3y 

T*  sitelz 

T  site2z 

T  siteZz 

(2.31) 


If  we  solve  for  p2  in  equation  (2.26)  we  can  see  that 


p2  —  A*2-\-  B2h2 

where 

A2  =  021^4  •  X  +  022^4  •  Y  +  CI23A  ■  Z ; 

B2  =  a,2iB  ■  X  +  CI22B  ■  Y  +  CI23B  ■  Z 

Similarly  from  (2.26)  if  we  solve  for  p\  and  p3  we  find, 

A*x  +  B*h2  Al  +  B*h2 

Pi  - - ;  P3  - - 

Cl  c3 

where 

Al,Bl,A^,B^  are  similar  in  form  to  A2  and  B2 

A\  =  —  ^cinA  •  X  +  cii2^4  ■  Y  +  ci\3A  ■  Z ^  ; 
B2  =  —  | a\\B  •  X  +  CL12B  ■  V  +  auB  ■  z'j 


(2.32) 


(2.33) 


(2.34) 


(2.35) 
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(2.36) 


Ag  —  —  •  X  +  Q.32^4  •  Y  +  033A  •  Z^j  ; 

5g  =  —  ^<231 B  ■  X  +  a32I?  ■  V  +  CZ33-B  ■ 

However  equations  (2.32)  and  (2.34)  still  contain  the  term  h2  defined  earlier  by  (2.17) 
which  is  still  unknown.  So  in  order  to  determine  the  magnitude  of  the  radius  vector 
to  our  target  at  the  central  time,  r2,  one  final  equation  is  needed. 


r2  —  rs«te2  +  {A  —  2*  +  i?2h2)2  +  2  (l2  •  rsite2^j  (A  —  2*  +  B2I12)  (2.37) 

Or,  in  its  more  recognizable  8th  order  polynomial  form, 


site  2 


+  A*2  +  2  (  L2  •  r  sitei]  A, 


—To 


2 A*BZn  +  2  L2  ■  rSJte2  B*ti  -  =  0 


(2.38) 


A  key  point  in  the  process  must  be  highlighted  here.  In  most  forms  of  the  Gauss 
algorithm,  the  emphasis  is  placed  on  selecting  the  largest  root  of  the  polynomial  as  the 
correct  value  of  the  magnitude  of  the  radius  vector  from  equation  (2.38).  This  is  true 
when  looking  out  from  the  origin  of  your  inertial  coordinate  frame  with  |r2|  >  \rsite2\. 
However,  in  this  problem  we  are  looking  in  towards  the  origin  |r2|  <  |rSjte2|  and  so 
root  selection  varies  slightly.  Selecting  the  largest  root  of  the  polynomial  will  return 
the  position  of  the  observing  sensor.  I  found  that  the  second  largest  positive  real  root 
is,  in  fact,  the  correct  value  for  the  magnitude  of  r2.  With  that,  it  is  then  possible  to 
use  this  method  to  find  three  position  vectors  for  the  target  satellite. 


2-4-2  Calculating  The  Velocity  Vector.  Using  the  Gauss  method  discussed 
in  the  previous  section,  it  is  possible  to  determine  the  position  vector  of  the  target 
satellite  at  three  different  times.  In  order  to  gain  a  complete  understanding  of  an  orbit 
in  three  dimensional  space,  the  velocity  vector  at  one  of  the  three  position  vectors 
must  also  be  known.  To  determine  the  satellite  velocity  vector  there  are  two  primary 
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techniques  which  are  used  depending  on  the  spread  of  data,  the  Gibbs  and  Herrick- 
Gibbs  methods.  The  Gibbs  solution  relics  on  the  geometry  of  the  problem  to  solve  for 
the  velocity  vector  and  works  best  when  the  spread  of  data  is  greater  than  five  degrees, 
whereas  the  Herrick-Gibbs  method  relies  on  a  Taylor  series  expansion  and  works  best 
when  the  spread  of  data  is  less  than  one  degree.  Between  one  and  five  degrees  both 
methods  work  reasonably  well  [15]  and  for  the  purposes  of  this  thesis  four  degrees  was 
selected  as  the  crossover  point  between  the  two  methods.  To  check  the  separation 
between  the  position  vectors  for  the  purposes  of  determining  the  appropriate  technique 
to  use  in  determining  the  velocity  vector  the  following  equations  were  used. 


COS  (CK12)  = 

r  1  •  r2 

(2.39) 

ri(  \f2\ 

cos  (a23)  = 

r2  ■  r3 

(2.40) 

\r2\  |r3| 

By  taking  the  sum  of  these  two  angles  we  can  determine  which  method  is  more 
appropriate  to  the  data  to  use. 

2.4.3  Gibbs  Method.  As  mentioned  earlier,  the  Gibbs  method  relies  on  a 
geometric  solution  to  solve  for  the  velocity  vector.  The  solution  method  presupposes 
that  we  know  three,  non-zero,  coplanar  position  vectors,  which  represent  three  se¬ 
quential  times  of  the  target  satellite  in  its  orbit.  From  the  Gauss  method  we  have 
the  appropriate  data  on  the  satellites  position  and  the  associated  sequential  times. 
When  using  this  method,  it  is  important  to  know  the  order  of  the  position  vectors  in 
time  because  it  relies  on  several  cross  products  in  its  solution.  We  begin  by  defining 
a  vector,  D,  as: 


D 


fi  x  r  2  +  r  2  x  r  3  +  r  3  x  -ri 


(2.41) 
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So  D  is  a  vector  which  is  perpendicular  to  the  plane  formed  by  the  three  position 
vectors  of  the  target  satellite.  We  also  need  to  define  a  second  vector  perpendicular 
to  the  plane,  N,  as: 

N  —  n  (r2  x  r3)  +  r2  (r3  x  fi)  +  r3  (fi  x  r2)  =  pD  (2.42) 

If  the  position  vectors  of  the  target  satellite  are  not  in  the  same  plane  then  the 
Gibbs  method  will  fail.  However,  for  non-maneuvering  spacecraft,  it  is  a  reasonable 
assumption  that  the  motion  is  planar.  Observational  errors  can  still  cause  the  vectors 
to  appear  non-planar.  A  method  of  solving  for  the  planar  separation  between  the 
position  vectors  is: 

Oicoplanar  =  90°  -  COS"1  (  ^  (2-43) 

V  r  2  X  r3|  |n|  J 

Since  the  observational  data  collected  from  the  focal  plane  array  contains  error,  this 
value  will  never  be  exactly  zero  so  it  is  necessary  to  choose  a  tolerance  for  this  value. 
Given  the  quality  of  the  data  analyzed,  a  tolerance  of  less  than  one  degree  proved  to 
be  sufficient.  Continuing  on  with  the  Gibbs  method,  we  need  to  define  one  additional 
vector  quantity  before  we  can  find  the  velocity  vector  at  the  central  time,  S. 

S  =  n  (r2  -  f i)  +  r2  (f3  -  f*i)  +  r3  (fq  -  r2)  (2.44) 

Finally,  we  can  calculate  the  velocity  vector  at  the  central  time  and  obtain  a  complete 
picture  of  the  orbit  at  our  selected  epoch  time. 

v~2  =  v2\liD0xi,2+]fii>s  (2’45) 

2-4-4  Herrick-Gibbs  Method.  The  second  method  which  was  used  to  solve 
for  the  velocity  vector  of  the  target  satellite  was  the  Herrick-Gibbs  method.  It  is 
used  when  the  vectors  are  too  close  to  one  another  for  the  Gibbs  method  to  work 
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successfully.  This  method  also  uses  the  three  observation  times  associated  with  the 
three  sequential  position  vectors  used  in  the  Gibbs  method.  As  in  the  Gibbs  method, 
it  is  necessary  to  conduct  a  coplanar  check  of  the  position  vectors  using  equation 
(2.43)  to  determine  the  validity  of  the  data  for  use.  Before  calculating  the  velocity 
the  change  in  time  between  position  vectors  of  the  target  can  be  defined  as: 

At3i  —  G  —  t] 

At32  —  t3  —  t2  (2.46) 

A^21  =  t2  —  t\ 

Now  using  a  truncated  Taylor  series  expansion  we  can  solve  for  the  velocity  vector  at 
the  central  epoch. 

Vo  =  — At39  (  — — - - 1 — \  G  +  (A t32  —  At2i)  ( -7 — 7 - 1 — \  G 

3“Va2iA  t31  12rf  /  1  ^J\A21At32  12  r\) 

(2.47) 

+At21  + Gi) ?3 

The  final  result  of  the  Gauss  method  is  an  estimate  of  the  position  and  velocity 
vectors  of  the  target  satellite.  These  values  must  be  converted  to  the  classical  orbital 
elements  used  in  the  state  vector  during  the  non-linear  least  squares  process  that  will 
be  discussed  in  Chapter  3.  All  of  the  classical  orbital  elements  are  calculated  from 
the  estimate  provided  by  Gauss,  but  only  the  four  values  in  the  state  vector  are  kept. 
The  full  element  set  calculated  by  Gauss  only  becomes  useful  when  the  target  is  not 
identified  in  the  AFSPC  catalog  as  discussed  in  Chapter  4. 

2.4.5  Spherical  Trigonometry.  The  second  initial  orbit  determination  tech¬ 
nique  utilizes  a  simpler  approach  to  obtain  an  incomplete  picture  of  the  targets  orbital 
elements,  but  it  is  enough  to  identify  it  within  the  catalog.  This  method  assumes  that 
the  observed  satellite  is  at  or  near  the  Earth’s  surface  in  a  circular  orbit.  This  as¬ 
sumption  does  introduce  error,  but  it  does  not  significantly  affect  the  information 
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Figure  2.5:  Initial  Orbit  Determination  Using  Spherical  Trigonometry 

about  the  orbital  plane,  which  is  defined  by  right  ascension  of  the  ascending  node 
(RAAN),  D  ,  and  inclination.  It  is  these  two  elements  that  are  initially  used  to  find  a 
match  in  the  AFSPC  catalog.  Using  some  simple  spherical  trigonometry  it  is  possi¬ 
ble  to  obtain  enough  information  to  begin  the  least  squares  iteration,  see  Figure  2.5. 
This  method  relics  on  knowing  the  location  of  the  reflected  signature  from  a  target 
relative  to  the  background  of  the  Earth.  Typically  this  is  how  fastwalkers  are  seen 
on  imagery  taken  from  a  platform.  After  being  processed  the  imagery  is  fitted  with  a 
latitude  and  longitude  grid.  Using  the  known  location  of  the  target  fastwalker  at  two 
closely  related  times  we  can  find  the  angle  between  true  North  and  the  orbit  plane, 
the  azimuth,  by: 


(2.48) 


It  is  then  possible  to  find  the  angle  A  in  Figure  2.5  if  we  define  latitude  at  epoch  as 
S  through, 


A  =  atari  ( sin  (<5)  tan  ( Az )) 


(2.49) 
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We  can  use  A  to  find  the  RAAN.  If  we  define  as  the  longitude  East  of  the  meridian 
and  Ogmst  as  the  Greenwich  mean  sidereal  time,  then  the  the  node  can  be  found  using 
equation(2.50). 

fl  =  Xe  +  Oqmst  —  A  (2.50) 

The  argument  of  latitude,  0  ,  can  then  be  found  using  the  latitude. 

6  =  acos  (cos  (A)  cos  (5))  (2.51) 


The  inclination  of  the  orbit  is  then  found  easily. 


sin  (5) 
sin  ( 0 ) 


(2.52) 


If  we  assume  a  small  change  in  time  between  observations  we  can  also  find  0  by, 


0  = 


\J (ALongitude)2  cos2  ( Lat )  +  (A Latitude)' 


to  —  1 1 


(2.53) 


Finally  using  0  we  can  obtain  an  estimate  of  the  satellite’s  altitude  above  the  Earth’s 
surface.  If  we  define  /i  as  the  gravitational  constant  of  Earth,  R,Earth  as  the  radius  of 
the  Earth,  and  h  as  the  altitude  we  can  find  altitude  using, 


0 


t* 


{.■O^Earth  T  h) 


(2.54) 


which  is  related  to  the  TLE  set  value  of  revolutions  per  day  by  Kepler’s  third  law. 
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III.  Non-Linear  Least  Squares 


The  method  used  to  refine  the  initial  state  vector,  found  through  either  the 
Gauss  initial  orbit  determination  or  using  the  spherical  trigonometry  technique,  is 
non-linear  least  squares.  The  initial  orbit  determination  techniques  discussed  earlier 
are  deterministic  in  their  approach.  Both  methods  use  an  equivalent  amount  of  known 
values  to  solve  for  a  set  of  equivalent  unknowns.  What  both  of  them  lack  is  an  ability 
to  account  for  less  than  perfect  data  and  for  additional  observations  outside  those 
necessary  to  obtain  the  unknowns.  Non-linear  least  squares  allows  for  the  introduction 
of  noise  into  the  data  as  well  as  for  the  incorporation  of  additional  observations  to 
compensate  for  that  noise. 

Least  squares  dates  back  to  the  work  of  Legendre  and  Gauss  at  the  beginning  of 
the  19th  century  [14],  Legendre  first  published  a  method  recognizing  the  importance 
of  observational  errors  and  proposed  the  least  squares  solution  as  a  method  for  min¬ 
imizing  those  errors  in  1806.  Gauss  later  published  a  probabilistic  method  for  using 
least  squares  in  1809,  but  claimed  to  have  developed  and  used  it  as  early  as  1795 
and  so  is  given  the  credit  for  having  originated  the  method.  Regardless  of  whoever 
originated  the  idea,  it  has  undergone  much  refinement  over  the  two  centuries  since  it 
was  first  developed. 

3. 1  General  Method  of  Non-Linear  Least  Squares 

Non-linear  least  squares  was  developed  by  Gauss  to  model  the  non-linearity  of 
real  world  systems.  For  the  purposes  of  this  thesis  we  have  a  non-linear  dynamic 
system  and  non-linear  observation  geometry  which  requires  the  use  of  this  method. 

Since  the  system  dynamics  are  available  as  an  explicit  solution 

x(t)  =  h(x(tQ)  ,t)  (3.1) 
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which  is  a  function  of  time  and  initial  conditions,  we  can  linearize  the  dynamics  in 
the  following  manner. 

Sx  (t)  =  <f>  (Mo)  Sx  (to)  (3.2) 

Where  the  quantity  Sx  will  be  small.  Before  the  dynamics  can  be  successfully  lin¬ 
earized  an  initial  trajectory  must  be  obtained  to  linearize  about.  The  ’true’  value  of 
the  system  state,  x0,  is  unobtainable  so  we  use  a  reference  trajectory  xref,  which  is 
hopefully  close  to  the  true  system  state.  The  reference  trajectory  is  obtained  from  us¬ 
ing  either  the  Gauss  initial  orbit  determination  technique  or  by  utilizing  the  spherical 
trigonometry  technique. 

The  condition  of  the  available  observational  data  must  also  be  considered.  This 
data  is  a  non-linear  function  of  the  system  state  at  the  current  time  and  the  observa¬ 
tion  geometry  and  can  be  written  as, 


Zi{ti)  =  G(x(ti),ti)  (3.3) 

The  observation  relation,  described  by  equation  (3.3),  can  be  linearized  in  a  manner 
similar  to  equation  (3.2).  Since  our  sensor  obtains  very  accurate  observations,  but 
not  perfect  observations  of  the  true  state,  it  does  not  produce  perfect  data,  z0.  So  we 
can  expect  that  as  the  error  in  the  observations  goes  to  zero  the  error  in  the  state  also 
approaches  zero.  If  the  error  is  small  enough  than  we  can  linearize  the  difference.  If 
we  first  describe  the  difference  between  the  true  state  and  calculated  state  as, 

x  =  Xq  +  Sx  (3.4) 

we  can  then  describe  the  difference  between  the  true  error  and  the  actual  data. 


e  =  z  —  Zq  =  G  (x,t)  —  G  (x0,t) 
=  G  (xq  +  5x,t)  —  G  (xo,t) 

“  aiSx{t) 


(3.5) 
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Since  the  observations  are  not  perfect  we  can  equate  the  true  error  to  the  residual  in 
the  data,  r.  The  residual  has  two  parts.  The  first  part  comes  from  the  error  in  the 
reference  trajectory  from  our  initial  orbit  determination  method  and  the  second  part 
from  the  errors  in  the  observations.  Using  our  reference  trajectory  we  can  compute 
the  observations  that  would  result  from  it. 


n  =  Zi-G  (xref  (■ ti ) ,  U)  (3.6) 

The  residual  vector,  r,,  is  a  measure  of  the  difference  between  the  observation  and 
the  expected  observation.  It  is  related  to  the  change  in  the  reference  trajectory,  5x, 
which  corrects  the  reference  trajectory  to  a  value  close  to  the  true  trajectory  of  the 
observed  target.  If  we  first  define  a  matrix,  Hi ,  which  relates  the  error  in  the  state 
observation  to  the  error  in  the  reference  trajectory  as 

dG 

Hi  (U)  =  —  (xref  (U) ,  U)  (3.7) 

The  residuals  in  the  observations  can  be  related  to  the  correction,  Sx  by 


r i  ~  H,i5x  (U)  =  Hi$  (t,  t0)  Sx  (t0) 


=  TiSx  (t0) 


(3.8) 


Finally,  a  covariance  matrix,  Q,  which  compensates  for  unit  compatibility  and  unit 
conversions  must  be  defined.  Assuming  each  measurement  taken  by  the  sensor  is 
independent,  Q  can  be  written  as 


/ 


Q  = 


ok 


(To 


0 

0 


0  0 


°N  ) 


(3.9) 
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where  a  is  the  statistical  variance  in  each  observation.  Combining,  we  can  write  that 
the  correction  to  the  reference  trajectory  can  be  written  as 


Sx(ta)  =  (TtQ-'T)  'r'Q  V 


(3.10) 


The  associated  covariance  can  then  be  written. 


P,r.  =  (TtQ-1T)  1 


(3.11) 


Using  equation  (3.10),  we  can  then  obtain  an  estimate  of  the  true  trajectory. 


X  (to)  =  Xref  {to)  +  Sx  (f0) 


(3.12) 


A  simple  modifcation  of  the  algorithm  in  Wiesel  [17]  can  best  describe  the  entire 
method. 

1.  At  each  observation  time  U,  propagate  the  calculated  state  vector  at  that  time 
to  the  epoch  time  using  the  state  transition  matrix,  <f>  (U,f0)- 

2.  Calculate  Ht  for  the  current  observation  as  well  as  the  residual  r*  =  Zi  —  G  {x). 

3.  Calculate  H j<b. 

4.  Add  this  T*  to  previous  ones  to  build  the  matrix  Tj Q~xTi  and  the  vector 
TJ Q- 1  r\  which  should  be  growing  larger  in  dimension  with  each  data  point. 
Repeat  these  steps  for  all  data  points  at  the  epoch  time. 

5.  Once  complete  with  all  data  points  calculate  the  covariance  P$x  by  using  equa¬ 
tion  (3.11). 

6.  Then  calculate  the  correction  to  the  reference  trajectory  using  equation  (3.10). 

7.  Correct  the  reference  trajectory  using  equation  (3.12). 
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8.  Determine  convergence  after  the  second  iteration  based  on  the  size  of  the  cor¬ 
rection.  ff  the  process  has  not  converged  continue  to  repeat  the  process  with 
the  newly  corrected  reference  trajectory  until  it  does  so. 

Each  of  the  several  functions  described  here  will  now  be  developed  in  detail. 

3.2  The  State  Transition  Matrix 

In  order  to  improve  upon  the  deterministic  reference  orbit  found  through  the 
initial  orbit  determination  techniques,  we  must  utilize  a  larger  quantity  of  data.  Least 
squares  allows  us  to  do  this  but  to  do  so  our  observations  must  be  propagated  to  a 
single  epoch  time.  Let  us  begin  by  defining  the  state  vector 

(3.13) 

Where,  as  before,  fl  is  the  right  ascension  of  the  ascending  node,  i  is  the  inclination 
of  the  orbit  and  9  is  the  argument  of  latitude,  which  is  equal  to  the  mean  anomaly  in 
a  perfectly  circular  orbit.  The  value  9  can  still  be  used  to  find  an  approximation  of 
altitude.  Since  all  the  values  are  constants  except  9,  which  is  a  function  of  time,  we 
can  define 

(3.14) 


X  =  F(X)  = 


0 

0 

9 


X  = 


n 

i 

9 

9 


0 


The  state  transition  matrix  is  then  simply, 


10  0  0 


(3.15) 


0  0  0  1 


3.3  The  Observation  Function,  G 

This  thesis  does  not  use  actual  platform  observations  so  they  must  be  generated 
separately.  The  observations  are  azimuth  and  elevation  as  measured  from  the  sensor 
platform  in  the  sensor  topocentric  frame.  Using  a  orbit  modeling  software  called 
Satellite  Tool  Kit  developed  by  Analytical  Graphics  Incorporated  we  can  obtain  the 
position  of  the  target  r  and  the  sensor  R  in  the  Earth  Centered  Inertial  frame.  The 
range  between  them  can  be  found  with 


(3.16) 


Pec  i  =  r  —  R 


This  range  then  needs  to  be  rotated  into  the  sensor  frame  in  order  to  simulate  our 
observations.  The  rotation  matrix  is  given  by 


(3.17) 


where  we  define 


(3.18) 


|fc  x  R | 


(3.19) 


S  =  Ex  Z 


(3.20) 
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and  k  =  0  0  1  is  the  unit  vector  in  the  Z  direction  in  the  ECI  frame.  The  range 

in  the  SEZ  frame  can  then  be  found  by 


Psez  —  Reci^sezPeci 


(3.21) 


Then  if  we  separate  out  the  components  of  the  range  in  the  SEZ  frame  as 


Ps 


Psez  —  pE 


(3.22) 


Pz 


we  can  write  our  observations  of  azimuth  and  elevation  in  terms  of  what  we  wish  to 
find,  the  range  between  our  target  and  sensor. 


(3.23) 


Az  =  atari 


(3.24) 


3-4  Measurement  Jacobi  Matrix,  H 

From  the  derivation  of  the  general  method  of  non-linear  least  squares,  in  par¬ 
ticular  equation  (3.7),  we  can  see  the  need  for  a  relationship  between  the  error  in  the 
state  and  the  error  in  the  reference  trajectory.  We  can  define  this  relationship  with 
regard  to  the  observations  (obs)  as 


dobs  _  dobs  dpSEZ  dp  eci 
dr  eci  dpsEZ  dpEci  dr  Eci 


(3.25) 


where 


(3.26) 
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which  is  the  rotation  matrix  we  found  previously  in  equation  (3.17)  and 

dpEci 


dr 


=  I 


ECI 


(3.27) 


is  simply  equal  to  an  identity  matrix.  The  challenge  here  is  the  first  term  which 
cannot  be  directly  related  but  must  be  further  broken  into, 


dobs  dobs  drECi 


dpsEZ  Steci  dX 


(3.28) 


In  order  to  find  the  relationship  between  the  position  vector  in  the  ECI  frame  and 
the  state,  which  we  will  call,  ,  we  first  need  to  define  a  rotation  between  the  ECI 
frame  and  the  orbital  frame  in  terms  of  the  state  vector.  If  we  first  define 


T pqw 


recos  ( 6 ) 
resin  ( 9 ) 

0 


(3.29) 


where  with  p  being  the  gravitational  constant  of  Earth  and  re  is  the  calculated  distance 
of  the  target  from  its  center. 
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We  can  than  find  teci  (W)  by  rotating  about  the  node  and  inclination. 

tec i  (X)  =  R3  (-12)  R3  (- i )  rp 

In  explicit  form 
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the  matrix  $  is  then 
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Finally,  we  need  the  function  ,  which  we  will  define  as  H  ■  Reci^sez- 
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H3 


We  can  then  find  and  build  T. 
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This  completes  the  functional  derivation  of  the  values  needed  to  execute  the 
non-linear  least  squares  algorithm.  The  actual  implementation  is  more  complicated 
and  as  a  reference  the  code  is  included  in  Appendix  2. 
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IV.  Data  and  Results 


4-1  Data  Accuracy  and  Precision 

An  important  factor  to  consider  in  the  development  of  any  model  for  a  real 
world  sensor  is  the  accuracy  of  its  measurements.  In  order  to  investigate  the  affect 
of  error  in  the  measuring  sensor  on  the  ability  of  the  method  being  developed  to 
accurately  identify  the  target  satellite,  a  range  of  reasonable  values  for  this  error 
must  be  determined.  The  primary  source  of  noise  in  the  data  was  assumed  to  be  the 
resolution  of  the  focal  plane  array  of  the  sensor  satellite,  known  as  ground  sample 
distance  (GSD). 

To  determine  the  effects  of  GSD  on  the  data  the  position  of  the  sensor  is  assumed 
to  be  known  to  a  level  of  precision,  which  for  our  purposes,  we  can  assume  to  be 
perfectly  known.  Another  possible  source  of  error  that  is  assumed  to  be  zero  is  known 
as  jitter,  this  refers  to  the  pointing  accuracy  in  the  sensor.  It  exists  in  any  man¬ 
made  optical  platform  but  is  small  compared  to  the  error  in  the  known  position  of  the 
sensor,  which  is  considered  to  be  perfectly  known.  Since  additional  sources  of  error  are 
smaller  or  have  no  effect  on  the  method  or  collection  of  data,  all  of  these  error  sources 
including  dead  pixels  and  dark  current  are  ignored.  Finally,  any  error  generated 
by  sub  pixel  processing,  centroiding  or  any  other  software  method  of  increasing  the 
resolution  of  the  focal  plane  array  is  considered  to  be  small,  and  other  than  changing 
the  value  of  the  GSD,  is  ignored. 

Any  satellite  in  a  highly  elliptical  orbit  spends  most  of  its  orbital  period  near  the 
apogee  of  its  orbit.  We  will  therefore  calculate  the  range  of  error  in  the  measurements 
of  azimuth  and  elevation  from  an  altitude  at  or  near  apogee  for  a  typical  HEO  orbit, 
36,000  kilometers.  The  range  of  reasonable  ground  sample  distances  that  will  be 
considered  for  investigation  is  one  kilometer  to  one  meter.  Ground  sample  distance 
is  simply  the  area  on  the  ground  that  a  sensor  sees  from  orbit  for  each  pixel  on  the 
focal  plane  array.  So  a  ground  sample  distance  of  one  meter  means  that  each  pixel  on 
the  focal  plane  array  of  the  sensor  sees  a  one  meter  by  one  meter  square  of  ground. 
While  many  civilian  sensors  working  in  the  visible  light  wavelengths  regularly  achieve 
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ground  sample  distances  of  one  meter  or  better,  the  technology  of  systems  operating 
outside  the  visible  wavelengths  lags  far  behind.  Calculating  the  range  of  this  error  is 
found  simply  by 


(4.1) 


Using  this  calculation  the  magnitude  of  the  errors  induced  by  GSD  is  seen  to  be  quite 
small  from  3.183e-3  to  3e-6  degrees.  This  level  of  error  is  unrealistically  small  for  real 
time  data.  Post  processing  of  data  can  significantly  improve  the  quality  of  the  data  to 
an  unknown  degree,  but  may  be  unnecessary  since  the  methods  developed  are  valid 
with  as  much  as  le-2  degrees  of  error  in  the  data. 

Jf.,2  Test  Case  Selection 

The  selection  of  test  cases  and  the  validation  of  the  method  was  done  by  ana¬ 
lyzing  the  satellite  catalog  for  inclination  and  right  ascension  of  the  ascending  node 
(RAAN) .  The  third  value  being  effectively  solved  for  in  the  least  squares  process  was 
the  argument  of  latitude.  This  value  is  a  measure  of  where  the  satellite  is  in  its  orbit 
which  is  constantly  changing  and  so  is  not  used  in  an  initial  search  of  the  catalog  for 
the  target  satellite.  To  use  the  argument  of  latitude,  each  item  in  the  catalog  would 
need  to  be  propagated  to  the  epoch  time  being  used.  This  is  simply  not  a  realistic 
method  of  conducting  an  initial  search  for  the  target  satellite,  but  can  yield  useful 
results  once  the  field  of  candidate  targets  is  narrowed  significantly.  Five  test  cases 
were  selected  using  satellites  from  various  bands  of  inclination.  Using  Figure  4.2  we 
can  clearly  see  bands  of  crowded  inclination  and  begin  a  selection  of  useful  test  cases. 
This  plot  represents  the  entire  catalog  of  satellites  as  generated  by  Air  Force  Space 
Command  as  of  28  September  2007  and  includes  11,231  objects  in  orbit  ranging  from 
fully  functional  satellites  to  small  pieces  of  debris. 

For  the  sensor  platform  in  HEO  orbit,  objects  below  40  degrees  of  inclination 
do  not  enter  the  useful  field  of  view  to  any  significant  degree,  so  we  will  focus  on 
objects  above  this  inclination  for  this  sensor  orbit.  The  three  largest  concentrations 
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of  satellites  can  be  seen  at  83,  74  and  98  degrees  of  inclination,  so  a  test  case  within 
each  of  these  bands  was  selected.  Finally,  one  test  case  was  selected  within  the  60-70 
degree  inclination  region  and  one  in  the  40-60  degree  region.  The  two  line  element 
sets  of  the  selected  test  cases  for  the  HEO  sensor  are: 

Test  Case  1  HEO:  98  degrees  inclination  LEO 

THOR  BURNER  2  R/B 

1  03271U  68042B  07269.97501982  .00000015  00000-0  29790-4  0  2392 

2  03271  098.8574  065.9071  0051791  293.5969  065.9783  14.15994498  29159 

Test  Case  2  HEO:  82  degrees  inclination  LEO 

COSMOS  2416 

1  28908U  05048A  07270.55302940  .00000017  00000-0  40682-4  0  5277 

2  28908  082.4733  359.0931  0005782  305.7944  054.2602  12.55368842  80898 

Test  Case  3  HEO:  74  degrees  inclination  LEO 

SL-8  R/B 

1  22676U  93036B  07270.31673174  .00000030  00000-0  20391-4  0  2997 

2  22676  074.0421  190.2875  0018551  003.8531  356.2753  14.33042505746561 

Test  Case  4  HEO:  66  degrees  inclination  LEO 

JASON 

1  26997U  01055A  07269.49083870  -.00000038  +00000-0  +10000-3  0  08486 

2  26997  066.0446  298.5464  0008096  264.5280  095.4809  12.80929963271384 

Test  Case  5  HEO:  51  degrees  inclination  LEO 

ISS  (ZARYA) 

1  25544U  98067A  07270.64846184  .00007688  00000-0  54207-4  0  607 

2  25544  051.6348  295.2535  0002708  045.3662  054.4413  15.75238918506947 

When  selecting  test  cases  for  use  with  a  sensor  in  GEO  orbit  there  are  several 
different  factors  that  need  to  be  investigated.  Unlike  a  HEO  orbit,  a  GEO  sensor  is 
much  more  likely  to  observe  a  target  that  is  also  in  GEO  orbit  and  very  close  to  the 
sensor,  so  this  case  must  be  investigated.  Additionally,  a  GEO  sensor  may  see  targets 
crossing  the  equatorial  plane  that  are  in  HEO  or  GEO  and  so  they  will  not  be  near 
the  surface  of  the  Earth  when  they  are  observed.  A  GEO  sensor  may  also  see  targets 
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in  highly  eccentric  orbits  and  the  ability  of  the  model  to  determine  their  orbits  must 
be  tested.  Inclined  LEO  orbits  and  polar  orbits  must  also  be  tested  to  show  that  the 
method  continues  to  work  with  these  orbits  as  it  does  when  the  sensor  is  in  a  HEO 
orbit.  Finally,  one  unique  orbit  is  tested,  Chinasat  7  is  a  satellite  with  a  28  degree 
inclination  that  has  a  near  GEO  altitude  when  it  crosses  the  equatorial  plane  and  its 
solution  is  investigated. 

Test  Case  1  GEO:  GEO,  0  degree  inclination  GEO 
SPACEWAY  3 

1  32018U  07036A  07268.23359173  -.00000182  +00000-0  +00000-0  0  00393 

2  32018  000.1973  288.5111  1772164  041.4316  020.7444  01.00364884000543 

Test  Case  2  GEO:  near  GEO,  28  degrees  inclination,  eccentric  orbit 
CHINASAT  7 

1  24282U  96048A  07264.78888470  -.00000229  00000-0  10000-3  0  1129 

2  24282  028.1690  345.6951  3099610  354.9022  002.6532  01.06657574  45729 

Test  Case  3  GEO:  high  altitude  apogee,  28  degrees  inclination,  eccentric  orbit 
N-2  R/B(2) 

1  12787U  81012C  07270.74053674  .00006018  00000-0  64665-3  0  1367 

2  12787  028.2407  310.6441  6437703  242.5277  039.3189  03.43267090274783 

Test  Case  4  GEO:  51  degrees  inclination  LEO 
ISS  (ZARYA) 

1  25544U  98067A  07270.64846184  .00007688  00000-0  54207-4  0  607 

2  25544  051.6348  295.2535  0002708  045.3662  054.4413  15.75238918506947 

Test  Case  5  GEO:  99  degrees  inclination  LEO 
OPS  4467  A 

1  00812U  64031A  07270.22049292  -.00000005  00000-0  21485-4  0  510 

2  00812  099.8399  266.8336  0005731  143.8671  216.2893  14.23220549243754 

Test  Case  6  GEO:  20  degrees  inclination  LEO 
DELTA  2  R/B 

1  22015U  92039B  07270.49272632  .00000760  00000-0  22968-4  0  3947 

2  22015  020.7183  320.6482  0076224  136.6272  224.0191  14.98105611855144 
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4-3  Data  Generation 

The  data  for  the  analysis  of  each  of  these  test  cases  was  generated  using  version 
7  of  the  Satellite  Tool  Kit  (STK).  Following  the  selection  of  a  suitable  test  case,  as 
described  in  the  previous  section,  a  generic  sensor  platform  was  generated  either  in  a 
GEO  or  HEO  orbit  depending  on  the  test  case  being  run.  The  GEO  sensor  orbit  had 
an  altitude  of  36,000  km  and  exactly  0  degrees  of  inclination,  while  the  HEO  orbit 
had  an  apogee  of  36,000  km  and  an  inclination  of  63.4  degrees.  The  test  case  was 
then  propagated  forward  from  the  data  listed  in  the  AFSPC  catalog  entry  until  both 
the  sensor  and  the  target  were  within  a  reasonable  line  of  sight  of  each  other  with  the 
Earth  as  a  background.  The  data  then  generated  was  the  latitude  and  longitude  of 
the  suborbital  path  of  the  target,  the  position  of  the  target  in  the  ECI  frame,  and  the 
position  of  the  sensor  in  the  ECI  frame.  The  suborbital  path  data  generated  by  STK 
was  perfectly  nadir  so  care  has  been  taken  to  explain  the  weakness  of  the  spherical 
trigonometry  initial  orbit  determination  technique  which  relies  on  this  perfect  data. 
The  positions  of  the  sensor  and  the  target  were  used  to  generate  the  Azimuth  and 
Elevation  measurements  from  the  sensor  to  the  target  as  described  in  Chapter  3.  With 
this  data,  as  well  as  the  time  at  which  each  measurement  was  made,  the  analysis  was 
conducted. 

4-4  Debris 

Debris  in  orbit  is  a  well  known  fact.  Tracking  and  accounting  for  debris  has 
become  a  significant  field  of  research.  The  obvious  danger  that  these  objects  pose  to 
functional  satellites  in  orbit  is  of  great  concern.  A  fundamental  assumption  of  this 
thesis  is  that  observed  fastwalkers  are  in  low  Earth  orbit.  The  possibility  that  debris 
could  represent  a  fastwalker  must  be  explored. 

Only  one  version  of  the  Space  Command  catalog  was  used  for  this  thesis.  This 
issue  of  the  catalog  was  for  the  period  between  26  and  28  September,  2007  and 
contained  11,231  total  entries.  Of  these  entries  6,294  were  classified  explicitly  as 
debris,  or  about  56  percent.  Defining  what  constitutes  debris  is  difficult.  Space 
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Command  currently  tracks  all  debris  larger  than  10  cm  in  diameter  which  provides 
a  bottom  limit  on  the  size  of  debris.  However,  the  upper  limit  remains  undefined. 
Rocket  bodies  and  other  large  objects  are  not  identified  as  debris  so  we  can  assume 
that  those  items  classified  as  debris  are  items  that  represent  a  piece  of  what  was  once 
a  single  tracked  object.  One  clear  and  obvious  example  is  the  Fengyun  1C  satellite. 
On  January  17,  2007  China  tested  an  anti-satellite  missile  on  this  satellite.  The  test 
was  successful  and  the  satellite  broke  into  many  smaller  pieces.  Over  2,000  pieces  of 
debris  from  this  single  satellite  are  being  tracked  in  orbit.  If  none  of  these  pieces  are 
smaller  than  10  cm  in  diameter,  then  we  can  safely  assume  that  not  very  many  of 
the  pieces  of  this  one  satellite  remaining  in  orbit  are  too  much  larger  than  that.  The 
possibility  that  one  of  these  many  small  pieces  of  debris  could  represent  a  fastwalker 
observed  by  our  sensor  must  be  investigated. 

A  fastwalker  is  defined  as  a  bright  object  moving  across  the  focal  plane  array 
of  an  Earth  observing  sensor  at  a  velocity  too  great  to  represent  an  object  within 
the  atmosphere.  The  key  word  in  this  definition  is  the  term  bright.  It  means  an 
object  that  is  reflecting  enough  sunlight  in  the  direction  of  the  sensor  to  be  visible 
against  the  background  of  the  Earth.  If  we  assume  that  the  objects  remaining  from 
the  Fengyun  1C  satellite  are  larger  than  10cm  in  diameter  but  are  irregularly  shaped 
then  we  can  begin  to  build  a  picture  of  their  behavior  in  orbit.  Since  these  objects 
are  no  longer  stabilized  in  any  way  they  will  rotate  about  their  major  axis  or  align 
in  the  orientation  which  produces  the  least  drag  from  the  atmosphere.  Outside  the 
atmoshpere,  if  a  portion  of  these  objects  have  a  reflective  surface  they  will  in  fact 
reflect  sunlight  when  exposed  to  the  sun  for  an  amount  of  time  dependent  on  the 
period  of  the  object’s  rotation  about  it’s  major  axis.  This  period  of  time  will  be 
very  brief  and  the  reflection  will  be  in  an  irregular  direction  as  the  orbit  of  the  debris 
changes  so  does  it’s  position  in  relation  to  the  sun.  It  is  therefore  possible  for  these 
objects  to  reflect  visible  light  in  the  direction  a  sensor,  but  it  is  unlikely  that  they 
will  do  so  for  a  length  of  time  necessary  to  develop  a  data  track  in  the  absence  of 
drag.  If  the  orbit  is  not  at  a  high  altitude  significant  drag  is  present  and  the  object 
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will  align  itself  in  the  orientation  which  produces  the  least  drag.  This  will  create 
a  situation  where  observation  is  possible  but  the  small  size  and  irregular  reflecting 
surface  of  debris  will  make  this  unlikely. 

As  mentioned  previously,  the  catalog  has  a  large  number  of  items  classified  as 
debris  listed.  The  surprising  fact  is  that  the  majority  of  these  listings  can  be  accounted 
for  as  pieces  of  debris  from  just  a  handful  of  satellites  like  the  Fengyun  1C.  Many  of  the 
items  listed  as  debris  must  also  be  the  result  of  debris  created  during  the  staging  and 
boosting  phases  of  certain  classes  of  rockets.  However,  the  rocket  bodies  themselves 
are  not  classified  as  debris  and  are  not  eliminated  as  possible  fastwalkers.  The  reason 
for  this  rests  in  the  size  and  regular  shape  of  a  typical  rocket  body.  Since  most 
rocket  bodies  will  be  cylindrical  and  unstabilized  they  will  rotate  about  their  major 
axis.  During  a  single  rotation  the  cylindrical  shape  will  reflect  light  once  towards  any 
arbitrary  point.  For  this  reason,  observers  will  see  them  not  as  a  constant  reflection 
but  they  will  appear  to  wink.  Due  to  this  ability,  as  well  as  the  surface  area  available 
to  reflect  sunlight,  rocket  bodies  cannot  be  discounted  as  possible  targets  as  we  can 
debris  fields  like  that  of  the  Fengyun  1C  satellite. 

4-5  Test  Case  Results 

This  section  will  discuss  the  results  using  the  test  cases  previously  mentioned 
and  the  models  previously  developed.  In  developing  data,  several  assumptions  about 
the  length  of  time  in  which  the  sensor  would  be  able  to  view  the  target  had  to  be 
made.  From  apogee,  a  one  degree  field  of  view  for  the  sensor  allowed  it  to  see  a 
circular  area  with  a  diameter  of  over  1350  kilometers  on  the  Earth’s  surface.  This 
was  assumed  to  be  a  reasonable  field  of  view.  Hence,  the  maximum  amount  of  time  a 
target  would  be  visible  was  for  about  four  minutes  assuming  the  target  was  detected 
and  tracked  early  in  the  field  of  view.  The  time  between  observations  was  set  at  one 
second,  although  any  time  period  between  observations  is  acceptable  as  long  as  they 
lie  within  the  field  of  view.  A  minimum  of  two  observations  was  required  for  the 
spherical  trigonometry  initial  orbit  determination  technique  and  three  observations 
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for  the  Gauss  technique.  In  order  for  the  least  squares  process  to  have  any  significant 
value  we  must  have  more  observations  than  are  used  in  the  initial  orbit  determination 
technique.  Additionally,  the  Air  Force  Space  Command  catalog  used  was  assumed  to 
be  the  one  published  nearest  the  time  of  the  observed  fastwalker.  Since  inclination  is 
essentially  constant,  but  the  node  changes  slowly  due  to  the  oblateness  of  the  Earth, 
they  can  be  considered  constant  for  the  duration  of  the  data  track.  These  values  are 
also  treated  as  constants  when  searching  the  AFSPC  catalog  so  it  is  necessary  to  use 
the  catalog  published  closest  to  the  epoch  time  of  the  observed  fastwalker  so  that  the 
error  in  node  made  by  this  assumption  is  small.  It  is  not  necessary  to  propagate  each 
TLE  in  the  catalog  to  the  epoch  time  being  considered  to  successfully  identify  the 
target. 

4-5.1  Test  Case  1  HEO  -  Thor  Rocket  Body.  The  method  of  using  the  slope 
between  two  points  and  the  geometry  of  the  Earth’s  surface  resulted  in  an  estimate  of 
the  state  vector  which  was  very  good  at  nadir,  an  elevation  of  -90  degrees,  for  three  of 
the  four  values  in  the  state  vector  with  two  data  points  one  second  apart.  The  state 
vector  consisted  of  four  elements:  right  ascension  of  the  ascending  node,  inclination, 
argument  of  latitude,  and  altitude.  This  method  resulted  in  an  estimate  of  the  node 
as  68.12  degrees  and  inclination  of  100.03.  From  the  two  line  element  set  in  the 
previous  section,  we  see  that  these  values  are  within  two  degrees  of  the  values  needed. 
The  value  of  the  argument  of  latitude,  as  discussed  previously,  is  not  a  constant  and 
need  only  be  used  to  check  that  the  satellite  selected  as  the  target  from  the  catalog 
is  in  fact  the  correct  satellite  by  propagating  the  element  set  to  the  epoch  time  of 
the  solution.  Altitude  is  necessary  for  the  least  squares  process,  but  the  accuracy 
of  its  value  is  not  critical.  However,  the  trigonometry  technique  produces  extremely 
poor  estimates  of  the  altitude  which  are  too  poor  to  allow  the  convergence  of  least 
squares.  While  increasing  the  time  between  observations  improves  the  estimate  of  this 
value,  it  greatly  decreases  the  accuracy  of  the  other  three  values.  Additionally,  the 
improvement  in  the  estimate  of  the  altitude  still  remains  relatively  poor  in  comparison 
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to  the  accuracy  of  the  other  three  values.  The  estimate  is  usually  more  than  several 
hundred  kilometers  off  of  the  actual  altitude.  Relying  on  an  earlier  assumption,  that 
the  viewed  satellite  is  in  a  low  earth  orbit,  an  average  value  for  these  satellites  can  be 
selected  which  is  adequate  for  least  squares  to  function.  A  value  of  500  km  was  used 
which  is  a  9  of  about  0.0011  radians  per  second.  The  values  for  the  performance  of 
the  spherical  trigonometry  technique  being  discussed  here  are  from  nadir  observations. 
Nadir  observations  are  ideal  data  for  this  method,  but  are  unlikely  to  ever  be  obtained 
by  a  real  platform.  This  method  suffers  greatly  from  off  nadir  angles  and  its  usefulness 
will  be  discussed  along  with  the  least  squares  process  and  its  ability  to  compensate 
for  poor  observations. 


A  Trig 


RAAN  =  68.1209°  /  =  100.0316°  9  =  57.6837°  9  =  0.00001  rad/  sec 


-A^Trig  -A Truth 


RAAN  =  2.2138°  I  =  1.1746°  9  =  -.6857° 


The  Gauss  initial  orbit  determination  technique  uses  the  same  angular  data 
necessary  for  the  least  squares  process.  It  can  be  adapted  to  use  the  same  latitude 
and  longitude  of  the  backgroimd  data  as  the  previous  method,  but  is  no  more  accurate 
with  that  data  than  using  simple  trigonometry.  When  adapted  to  using  azimuth  and 
elevation  data,  the  Gauss  technique  eliminates  the  inability  to  directly  observe  the 
altitude  of  the  target  as  a  problem.  This  method  provided  an  excellent  estimate  of 
all  four  states.  The  estimated  values  of  the  state  vector  were: 


X, 


Gauss 


RAAN  =  67.6875°  I  =  99.0811°  9  =  50.8519°  9  =  0.00108  rad/ sec 


A  Gauss  A  Truth 


RAAN  =  1.7804°  I  =  .2241°  d  =  -6.1461c 


It  can  be  quickly  seen  from  the  correct  values  that  Gauss  has  given  us  an  estimate 
of  inclination  accurate  to  within  one  degree.  It  is  not  as  obvious  that  the  value 
of  altitude  is  also  extremely  accurate  in  comparison  to  the  previous  technique.  The 
difference  between  the  true  altitude  of  the  target  and  the  calculated  value  was  17.9  km. 
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Figure  4.2:  Azimuth  Residuals  With  3.18e-3  Noise  for  Thor  Rocket  Body 

In  comparison  to  traditional  techniques  of  satellite  tracking  using  radar  to  calculate 
range,  this  is  very  poor  but  for  our  purposes  it  is  more  than  adequate. 

The  purpose  of  the  non-linear  least  squares  process  is  to  account  for  additional 
observations  given  an  initial  estimate  of  the  state  in  order  to  reduce  the  effect  of  sensor 
errors.  The  accuracy  of  that  state  estimate  is  important  in  that  a  poor  estimate  of 
the  initial  state  vector  cannot  always  be  corrected  for  by  the  least  squares  process. 
On  the  other  extreme,  with  a  very  accurate  estimate,  it  does  not  always  improve  on 
the  initial  guess.  This  is  because  the  developed  non-linear  least  squares  relies  on  a 
two-body  formulation  of  the  system  dynamics,  which  was  judged  to  be  more  than 
adequate  for  the  short  time  period  over  which  observations  occur.  Using  the  results 
from  Gauss  as  our  initial  guess  the  least  squares  process  returns: 


Xnlls  — 


RAAN  =  65.3317°  I  =  99.1120°  9  =  52.0130°  9  =  0.000961  rad/sec 


Xnlls  ~  XTruth  — 


RAAN  =  -.7830°  /  =  .2130°  9  =  -4.985° 


These  values  were  obtained  with  random  noise  in  the  data  with  a  standard  deviation 
of  3.18e-3  degrees  which  was  shown  earlier  to  represent  a  one  kilometer  ground  sample 
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Elevation  Residuals 
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Figure  4.3:  Elevation  Residuals  With  3.18e-3  Noise  for  Thor  Rocket  Body 

distance.  The  difference  from  truth  represents  the  values  that  would  be  searched  for 
in  the  catalog,  meaning  that  the  -.7830  degree  difference  in  RAAN  represents  the 
difference  between  the  listed  value  and  that  in  the  catalog,  not  the  difference  between 
the  RAAN  at  epoch  and  this  value  which  is  smaller.  However,  the  difference  value  for 
6  does  represent  the  error  from  the  6  at  epoch  since  this  value  is  otherwise  meaningless. 
The  residuals  can  be  seen  to  be  very  small  in  Figure  4.2  and  Figure  4.3.  The  regular 
shape  of  the  plot  shows  that  model  error  is  much  more  significant  than  sensor  data 
error.  Also,  using  the  root  mean  square  as  the  convergence  criteria,  we  can  see  the 
rate  of  convergence  in  Figure  4.4  as  well  as  the  magnitude  of  the  corrections  being 
made  to  the  initial  state.  The  least  squares  process  converged  in  33  iterations  in  this 
scenario  and  as  noise  increased  this  number  went  up  slowly,  staying  below  100  for  a 
successful  convergence.  One  factor  to  note  is  that  the  good  initial  guess  resulted  in  a 
least  squares  data  run  which  was  only  slightly  more  accurate  than  the  initial  guess. 

The  covariance  matrix  is  very  small,  due  to  the  small  errors  associated  with 
the  observations.  R  was  found  that  the  covariance  was  not  a  good  predictor  of  the 
actual  error.  This  is  consistent  with  the  fact  that  the  least  squares  dynamics  model 
used  was  not  as  accurate  as  the  model  used  to  generate  data.  As  a  result,  the  search 
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Residual  Norm  vs  iterations 
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Figure  4.4:  Convergence 


of  the  catalog  for  objects  matching  the  estimate  had  to  look  well  beyond  the  range 
suggested  by  the  covariance  values. 
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Taking  the  values  from  the  spherical  trigonometry  initial  orbit  determination 
technique  that  were  more  than  two  degrees  off  with  perfect  nadir  data,  we  find  that 
the  least  squares  process  with  same  amount  of  noise  in  the  data  returns  the  same 
results  for  the  value  of  the  state  as  seen  above  with  the  Gauss  inputs.  Investigating 
the  data  sets  generated  for  this  thesis  for  a  sensor  in  HEO,  the  largest  off  nadir  angle 
of  elevation  was  found  to  be  5.2  degrees.  Since  the  sensor  cannot  look  at  an  angle 
too  far  off  nadir  before  creating  significant  orthorectification  effects  on  the  imagery 
from  the  curvature  of  the  Earth,  we  can  assume  a  10  degree  off  nadir  look  angle  as 
being  the  greatest  allowable.  Using  this  angle  the  largest  amount  of  error  possible  in 
both  latitude  and  longitude  is  less  than  one  degree,  which  will  translate  into  nearly 
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the  same  amount  of  error  in  the  values  of  node  and  inclination  for  the  target  satellite. 
The  least  squares  process  was  found  to  easily  compensate  for  as  much  as  5  degrees 
of  random  error  in  the  first  three  states,  when  the  estimate  of  altitude  remained  at  a 
fixed  500  km  and  the  level  of  noise  in  the  data  at  3.18e-3  degrees. 

The  results  when  we  lower  the  GSD  to  one  meter  are  worth  noting.  Both  Gauss 
and  the  least  squares  process  return  answers  that  are  more  accurate  than  with  a  one 
kilometer  GSD  but  these  improvements  are  small  for  the  least  squares  process.  Gauss 
improves  signihcnatly  when  used  with  more  accurate  data.  We  can  see  below  that  the 
results  are  nearly  similar  for  Gauss  and  least  squares.  What  becomes  more  apparent 
when  the  GSD  is  lowered  to  one  meter  is  shown  by  the  plots  of  the  residuals.  The 
regular  shape  of  the  residuals  shown  in  Figure  4.5  and  Figure  4.6  is  the  result  of  model 
error  being  greater  than  the  error  due  to  noise  in  the  data.  This  helps  explain  the 
poor  results  from  the  covariance  matrix  which  fails  to  account  for  the  model  error,  a 
significant  factor  that  is  not  readily  apparent  when  the  noise  in  the  data  is  3.18e-3 
degrees. 


X 


Gauss 


RAAN  =  66.1549°  /  =  98.8411°  6  =  51.4038°  0  =  0.00103  rad/ sec 


Xnlls 


XQauss  X Truth 


RAAN 


RAAN  =  65.8516°  /  = 


Xnlls  ~  Xxruth  — 


RAAN 


=  .2478°  /  =  -.0159° 
99.0382°  6  =  51.7056° 
=  -.0555°  I  =  .1808° 


6 

6 

e 


-5.5942° 
0.000999  rad/ sec 
-5.2924° 


Concerning  noise,  one  final  look  at  the  robustness  of  the  non-linear  least  squares 
process  is  warranted.  Due  to  requests  from  the  sponsor,  the  effect  of  ground  sample 
distance  was  investigated  with  all  other  sources  of  error  assumed  to  be  small  in  com¬ 
parison.  If  these  assumptions  are  relaxed  and  the  level  of  random  noise  is  raised  to 
10e-3  degrees,  representing  a  GSD  of  12.65  km,  the  least  squares  process  was  found  to 
still  converge  to  a  reasonable  guess  at  the  state,  allowing  for  a  match  to  be  found  in 
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Azimuth  Residuals 


Figure  4.5: 


Figure  4.6: 


Azimuth  Residuals  With  3e-6  Noise  for  Thor  Rocket  Body 


Elevation  Residuals  With  3e-6  Noise  for  Thor  Rocket  Body 
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Azimuth  Residuals 


Figure  4.7:  Azimuth  Residuals  With  10e-3  Noise  for  Thor  Rocket  Body 

the  catalog.  The  limiting  factor  at  this  level  of  noise  is  the  use  of  Guass  as  an  initial 
orbit  determination  technique,  10e-3  represents  the  error  point  at  which  Gauss  would 
begin  returning  unreasonable  estimates  of  the  slant  range. 


Xnlls  — 


RAAN  =  62.3955°  /  =  100.0189°  9  =  53.9412°  9  =  0.000749  rad/ sec 


XnLLS— Truth, 


RAAN  =  3.5116°  I  =  1.1615°  9  =  -3.0568° 


While  still  small,  the  residuals  show  again  that  the  level  of  noise  in  the  data  has 
overcome  the  model  error,  see  Figure  4.7  and  Figure  4.8. 

The  final  step  in  validating  the  search  for  the  target  is  to  actually  determine 
if  the  data  developed  by  the  methods  discussed  above  is  accurate  enough  to  identify 
the  target  satellite.  As  mentioned  above,  the  covariance  matrix  was  extremely  small 
so  we  need  to  begin  looking  in  a  much  larger  range.  If  we  look  simply  at  inclinations 
within  one  degree  of  our  estimated  value  from  our  results  with  3.18e-3  degrees  of  noise 
a  search  of  the  catalog  returns  2,276  matches.  If  we  include  the  right  ascension  of  the 
ascending  node  in  that  search  we  find  only  one  match,  which  is  our  target  satellite. 
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Elevation  Residuals 


Figure  4.8:  Elevation  Residuals  With  10e-3  Noise  for  Thor  Rocket  Body 

In  this  case  it  proved  unnecessary  to  propagate  the  orbit  to  find  that  this  is  indeed 
the  correct  target,  but  doing  so  confirms  that  it  is. 


4-5.2  Test  Case  2  HEO  -  Cosmos.  This  test  case  dealt  with  the  82  degree  in¬ 
clination  orbit.  Both  initial  orbit  determination  techniques  worked  fairly  well,  with  the 
Gauss  technique  once  again  proving  to  be  far  superior  than  the  spherical  trigonometry 
technique  even  when  provided  with  nadir  data.  The  spherical  trigonometry  technique 
continued  to  have  extremely  poor  estimates  of  the  target’s  altitude.  When  a  search 
of  the  catalog  was  conducted  three  objects  were  found  to  be  close  matches  in  both 
inclination  and  RAAN.  After  propagating,  the  Cosmos  satellite  was  easily  identified 
by  the  argument  of  latitude  at  epoch. 


A  Trig 


RAAN  =  361.839°  I  =  83.999°  9  =  125.469°  9  =  0.00000895  rad/ sec 


X, 


Gauss 


RAAN  =  358.2412°  /  = 


Xnlls  — 


RAAN  =  359.1751°  /  = 


Xnlls  ~  XTruth  — 


RAAN 


81.6152° 
82.1518° 
=  .2471° 


9  =  122.0505°  0  =  .000828  rad/ sec 
9  =  121.9621°  6  =  .000901  rad/ sec 
I  —  —.3162°  0  =  .O111° 
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4-5.3  Test  Case  3  HEO  -  SL-8  Rocket  Body.  This  test  case  was  in  the  74 
degree  inclination  band.  Here  least  squares  improved  on  Gauss  significantly  and  the 
result  was  very  close  to  the  actual  values.  The  search  of  the  catalog  resulted  in  three 
matches  within  one  degree  of  these  values.  Only  one  of  which  matched  after  including 
argument  of  latitude  following  propagation  to  epoch. 


X Trig 


RAAN  =  191.3505°  /  =  74.6582°  9  =  65.5391° 
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Gauss 


RAAN  =  189.3430°  I  - 


NLLS 
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Xnlls  —  X 


Truth  — 


RAAN 


74.3784° 
74.4999° 
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9  =  58.9821° 
9  =  58.9529° 
/  =  -.3311° 


9  =  0.00000897  rad/  sec 

9  =  0.00099  rad/ sec 
9  =  0.001005  rad/ sec 
9  =  -6.932° 


4-5-4  Test  Case  4  HEO  -  Jason.  The  inclination  tested  here  was  66  degrees. 
In  this  test  case  the  Gauss  technique  provided  a  poor  answer  in  both  RAAN  and 
inclination  but  the  least  squares  process  was  easily  able  to  converge  on  the  correct 
solution.  As  the  catalog  was  searched  outward  from  the  least  squares  results,  the  first 
match  proved  to  be  the  correct  target  after  propagating  to  confirm  the  argument  of 
latitude  at  epoch. 


Xj^rig 


RAAN  =  295.9209°  /  =  67.1962°  9  =  122.2495° 
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NLLS  ~ 
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XnLLS  ~  XTruth 


It  TTY 


62.7412° 
65.9347° 
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9  =  113.8477° 
9  =  115.1594° 
/  =  -.0643° 


9  =  0.00000790  rad/ sec 

9  =  0.001165  rad/sec 
9  =  .0009512  rad/ sec 
9  =  -6.193° 


4-5.5  Test  Case  5  HEO  -  International  Space  Station.  This  test  case  repre¬ 
sented  the  lowest  inclination  tested  at  51  degrees.  It  begins  to  show  a  slight  failing  of 
the  use  of  azimuth  and  elevation  data  as  the  portion  of  the  observed  orbit  is  nearly  flat 
with  respect  to  the  Equator  and  low  on  the  sensor  horizon.  This  resulted  in  only  small 
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changes  in  the  recorded  values  for  elevation,  which  resulted  in  the  poor  performance 
of  least  squares.  In  this  case  the  spherical  trigonometry  technique  proved  to  be  su¬ 
perior  in  RAAN  and  inclination  though  the  least  squares  and  Gauss  techniques  both 
still  worked  well  enough  to  identify  the  target.  A  search  of  the  catalog  around  these 
values  revealed  five  objects.  Two  of  these  objects  were  debris  from  the  Space  Station 
itself  and  could  be  ignored.  The  other  two  did  not  fit  when  propagated  forward.  The 
final  object  was  the  target. 


A  Trig 


RAAN  =  294.1198°  I  =  51.6595°  9  =  81.5910°  9  =  0.00000713  rad/ sec 
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RAAN  =  .2042°  / 


73.9015°  9  =  0.001132  rad/ sec 
73.7151°  9  =  .001125  rad/sec 
.3662°  9  =  -14.115° 


4-5.6  Test  Case  1  GEO  -  Spaceway  3.  When  using  a  sensor  in  GEO  the 
problem  of  identifying  a  target  becomes  more  difficult.  A  key  assumption  that  was 
made  when  working  with  a  sensor  in  HEO  no  longer  applies.  No  longer  can  we  assume 
that  the  object  being  viewed  is  in  a  LEO  orbit.  We  must  also  consider  whether  the 
object  is  also  in  GEO  or  in  a  Medium  Earth  Orbit  (MEO).  Spaceway  3  is  a  GEO 
satellite  that  allowed  us  to  explore  the  problem  of  another  satellite  passing  directly  in 
front  of  a  GEO  sensor  and  here  we  see  the  method  fail.  The  spherical  trigonometry 
technique  results  appear  to  be  close  in  inclination,  but  it  must  be  remembered  that 
those  results  were  generated  using  perfectly  nadir  data.  Nadir  data  is  extremely  un¬ 
likely  in  any  scenario  and  becomes  even  more  so  as  the  distance  between  the  observer 
and  the  target  decreases.  The  spherical  trigonometry  method  remained  valid  only 
when  the  target  was  near  the  Earth’s  surface  and  the  sensor  was  in  a  high  enough 
altitude  orbit  so  as  to  make  errors  from  viewing  at  non-nadir  angles  small.  When 
viewing  an  object  very  near  the  sensor,  these  non-nadir  angles  can  span  the  entire 
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region  of  space  visible  to  the  sensor. 


A Trig 


RAAN  =  -293.188°  I  =  0.1923°  9  =  70.8042°  9  =  0.00002796  rad/ sec 


A  Trig  At  ruth 


RAAN  =  -221.699°  I  =  0.0005°  6  =  1.541° 


Identifying  the  altitude  of  the  target  is  impossible  with  the  spherical  trigonom¬ 
etry  technique  but  Gauss  gives  us  an  estimate  of  the  target  satellite’s  altitude  that  is 
at  least  accurate  enough  to  determine  that  it  is  also  in  a  GEO  or  near  GEO  orbit.  In 
this  case,  using  Gauss,  we  can  determine  that  the  altitude  of  the  target  was  greater 
than  30,000  kilometers,  making  it  a  target  which  we  will  be  unable  to  identify. 


Gauss 


RAAN  =  182.2929°  /  =  171.2489°  6  =  171.5152°  0  =  0.00009442  rad/ sec 


X 


Gauss— Truth 


RAAN  =  -106.218°  /  =  8.5538°  9  =  99.1702° 


The  problem  of  attempting  to  identify  satellites  in  GEO  from  GEO  was  attempted  in 
much  greater  detail  in  Osedacz  [9].  One  difference  must  be  noted.  Osedacz  declared 
that  it  was  impossible  to  determine  if  the  object  being  viewed  was  in  GEO  or  LEO. 
The  Gauss  method  used  here  was  able  to  provide  an  estimate  of  the  position  vector 
of  the  target  that  was  more  than  adequate  to  determine  the  difference  between  a 
GEO  and  a  LEO  satellite.  Since  we  can  differentiate  the  two  using  Gauss,  either  of 
the  initial  orbit  determination  techniques  provides  an  adequate  enough  estimate  of 
the  state  vector  for  a  LEO  target  to  apply  the  least  squares  process.  In  fact,  while 
the  method  performs  poorly  with  all  satellites  in  GEO  orbit,  it  performs  very  well 
with  low  inclination  circular  orbits  near  LEO.  Finally,  when  attempting  to  use  the 
initial  estimates  of  the  state  calculated  above  for  the  target  in  GEO  the  least  squares 
algorithm  failed  to  converge.  Based  on  this  case,  it  is  clear  that  the  methods  being 
developed  in  this  thesis  cannot  be  used  to  identify  GEO  or  near  GEO  satellites. 
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4-5.7  Test  Case  2  GEO  -  Chinasat  7.  Chinasat  7  is  a  unique  satellite  in  that 
it  was  the  highest  inclined  orbit  with  a  GEO  altitude  in  the  catalog  that  was  used.  It 
was  selected  to  test  the  response  of  Gauss  and  the  least  squares  process  to  an  inclined 
GEO  orbit.  Once  again  the  spherical  trigonometry  technique  is  unreliable  for  the 
reasons  stated  earlier  but  Gauss  was  somewhat  surprising.  The  Gauss  method  failed 
to  identify  the  satellite  as  being  in  an  elliptical  orbit  and  found  the  satellite  to  have 
a  parabolic  path.  With  the  failure  of  both  initial  orbit  determination  methods  the 
use  of  the  least  squares  process  is  no  longer  possible.  The  Chinasat  7  test  case  shows 
that  identifying  any  object  too  close  to  the  sensor  with  angles-only  data  is  extremely 
difficult  no  matter  what  inclination  it  is  at. 

4-5.8  Test  Case  3  GEO  -  N-2  R/B(2).  This  test  case  represented  a  target 
satellite  with  a  relatively  low  inclination,  a  high  altitude  apogee,  and  large  eccentricity. 
At  the  point  in  which  it  crosses  the  equatorial  plane,  the  satellite  is  close  to  LEO 
altitudes,  but  still  several  thousand  kilometers  from  the  surface  of  the  Earth.  In  this 
test  case,  the  method  barely  succeeds.  As  you  can  see  below,  the  least  squares  process 
returned  an  answer  that  was  actually  farther  from  the  truth  than  Gauss  due  to  the 
failure  of  the  circular  orbit  assumption  upon  which  the  least  squares  method  relies. 
Using  Gauss,  we  can  quickly  identify  it  but  with  the  poor  estimate  of  inclination 
provided  by  the  least  squares  process  15  possible  targets  had  to  be  propagated  before 
the  target  satellite  could  be  identified  by  the  mean  anomaly. 


A  Trig 


RAAN  =  309.756°  7  =  30.6425°  0  =  175.8093°  9  =  0.000797  rad/ sec 
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RAAN  =  313.8226°  I 


RAAN  =  309.4777°  I 


NLLS  ~  XTruth 


RAAN 


27.7682° 
34.0038° 
=  -.6063° 


9  =  169.1596° 
9  =  176.0905° 
7  =  5.7978° 


9  =  0.0009025  rad/ sec 
9  =  0.0007839  rad/ sec 
9  =  1.0965° 
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4-5.9  Test  Case  4  GEO  -  ISS(Zarya).  In  this  test  case  Gauss  provided  us 
with  a  better  estimate  of  the  true  value  of  inclination  than  the  least  squares  process, 
but  the  target  was  easily  identified  using  either  method.  For  the  HEO  test  cases,  the 
ISS  represented  the  lowest  inclination  test  case  and  in  that  situation  we  saw  Gauss 
and  least  squares  perform  equally  well.  With  accurate  data,  Gauss  has  proved  to  be 
an  excellent  initial  orbit  determination  technique,  but  it  is  deterministic,  and  a  single 
piece  of  poor  data  would  greatly  affect  its  result. 
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RAAN  =  294.224°  I  =  54.2655°  6  =  19.1739°  9  =  0.001088  rad/ sec 
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RAAN  =  292.9583 
RAAN  =  292.8512 
Xnlls  ~  XTruth  = 


°  I  =  51.6694° 

°  /  =  52.0571° 
RAAN  =  -.0161 


d 

9 


=  14.1904° 
=  22.1802° 
/  =  .448° 


9 

9 

9 


0.001176  rad/ sec 
0.001314  rad/ sec 
2.188° 


4-5.10  Test  Case  5  GEO  -  OPS  4457  A.  In  this  test  case  we  see  significant 
improvement  using  the  least  squares  process.  This  follows  with  what  we  saw  in  the 
HEO  test  cases.  The  target  was  easily  identified  with  the  given  state  vector.  The 
ability  of  this  method  to  identify  the  target  has  been  aided  by  the  fact  that  the 
method  performs  very  well  at  higher  inclinations  where  the  orbits  are  much  more 
crowded  than  at  lower  inclinations  but  tend  to  be  more  circular. 


A  Trig 


RAAN  =  268.6538°  /  =  103.1968°  9  =  22.7718°  9  =  0.001044  rad/sec 
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100.0285°  9  =  15.0649°  9  =  0.001011  rad/sec 
98.6737°  9  =  20.4633°  9  =  0.000931  rad/sec 
.3515°  1  =  -1.119°  9  =  -1.6118° 
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4-5.11  Test  Case  6  GEO  -  DELTA  2  R/B.  This  test  case  shows  the  ability 
of  both  initial  orbit  determination  techniques  and  the  least  squares  process  to  work 
well  with  a  LEO  satellite  at  low  inclination.  Using  Gauss,  the  approximate  altitude 
of  the  orbit  can  be  used  to  determine  that  the  target  is  not  at  or  near  GEO.  Once 
it  has  been  determined  to  be  in  a  LEO  orbit,  the  target  can  be  identified.  There  are 
only  49  satellites  in  LEO  at  inclinations  below  28  degrees,  the  highest  inclination  at 
which  we  are  likely  to  see  an  object  near  GEO.  This  is  a  small  number  of  satellites 
in  comparison  with  the  total  number  of  objects  being  tracked  from  0  to  28  degrees  of 
inclination,  ft  is  likely  that  because  most  objects  in  this  inclination  rage  are  not  in 
LEO  most  observed  fastwalkers  in  this  inclination  range  will  be  unidentifiable. 


X 


Trig 


RAAN  =  317.111°  I  =  22.0569°  9  =  15.7482°  9  =  0.00101  rad/sec 


XQauSS 


Xnlls 


=  RAAN  -- 
=  RAAN  -- 
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316.3434° 

=  RAAN 


I  =  20.6494°  9  =  9.6695°  9  =  0.00114  rad/ sec 
I  =  20.6253°  9  =  15.9865°  9  =  0.0000103  rad/ sec 
=  -4.3046°  /  =  -.0927°  9  =  -1.4085° 


4-5.12  Direct  Orbit  Prediction.  It  has  been  made  clear  several  times  that 
the  general  method  being  developed  is  dependent  entirely  on  the  assumption  that  the 
target  object  is  located  in  the  catalog.  However,  since  mistakes  and  omissions  can 
be  made,  what  can  be  done  in  the  case  where  the  object  is  not  located  within  the 
catalog  to  identify  its  orbit?  We  have  developed  Gauss  to  provide  us  with  a  specific 
state  vector,  but  it  can  also  provide  us  with  an  estimate  of  the  orbital  elements  of 
the  target  object.  As  has  been  shown  in  the  test  cases  previously  mentioned,  Gauss 
can  provide  us  with  results  that  are  accurate  enough  to  determine  the  identity  of  the 
target  using  the  catalog  and  given  a  data  track  spread  over  more  than  30  seconds  for 
a  LEO  object,  about  .5  degrees  of  arc.  Gauss  however,  provides  us  with  a  poor  direct 
estimate  of  the  orbital  position  of  a  satellite. 
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ISS  Orbital  Element  Comparison 

Element 

ISS  Actual 

ISS  Gauss 

Difference 

Semi-major  Axis  (km) 

6716.786 

6760.078 

-43.292 

Eccentricity 

0.000482 

0.008314 

-.007832 

Inclination  (deg) 

51.574 

51.545 

.029 

RAAN  (deg) 

293.557 

293.530 

.027 

Arg  of  Perigee  (deg) 

87.830 

90.634 

2.804 

Mean  Anomaly  (deg) 

354.044 

343.568 

10.476 

Table  4.1:  ISS  Orbital  Element  Comparison 


Table  4.1  is  an  example  of  the  quality  of  the  orbital  elements  estimate  provided 
by  Gauss.  These  results  show  that  the  largest  error  lies  in  the  argument  of  perigee 
and  the  mean  anomaly.  The  net  effect  of  these  errors  is  that  the  predicted  position 
of  the  target  satellite  will  lag  far  behind  that  of  the  actual.  Initially,  there  are  930 
km  separating  the  calculated  versus  actual  position,  though  they  are  following  very 
similar  suborbital  paths.  This  error  grows  with  each  orbit  as  the  error  in  the  periods 
between  the  predicted  and  actual  satellite  orbits  increases.  After  one  full  day,  the 
error  between  the  predicted  and  actual  position  is  nearly  1/4  of  an  orbit.  While  the 
direct  calculation  of  the  orbit  does  not  provide  us  with  an  accurate  position  of  the 
target,  it  does  provide  enough  information  to  reacquire  the  target  after  only  a  few 
orbits.  As  the  entire  method  relies  on  the  fact  that  very  few  satellites  share  both  an 
inclination  and  a  RAAN,  it  would  be  possible  to  reacquire  the  target  and  develop 
additional  tracking  data  to  refine  the  orbit.  While  this  would  require  the  dedication 
of  satellite  resources  to  reacquiring  the  target,  it  would  be  possible  with  an  additional 
data  track  to  develop  a  far  superior  direct  calculation  of  the  orbital  elements  by  using 
either  the  methods  developed  by  Gooding  [5]  or  Osedacz  [9]. 
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V.  Conclusions 


Fastwalker  orbit  determination  can  be  performed  with  a  single  data  set  from  a 
HEO  orbit  and,  with  some  limitations,  a  GEO  orbit.  Contrary  to  the  broad  assertions 
in  a  previous  thesis  on  the  topic  that  orbit  determination  of  fastwalkers  is  not  possible 
with  a  single  data  set  [9],  it  is  possible  to  characterize  LEO  orbits  from  a  sensor 
platform  in  a  HEO  and  GEO  orbit  using  orbit  matching.  Data  collected  from  future 
occurrences  can  be  used  to  determine  the  orbit  of  the  target  object. 

Two  different  initial  orbit  determination  techniques  were  discussed.  Both  were 
successful  in  determining  an  estimate  of  the  state  vector  necessary  for  additional 
analysis.  The  Gauss  method  was  far  more  accurate  than  the  spherical  trigonometry 
technique  over  data  sets  that  spanned  several  minutes,  but  as  the  observation  time 
decreased  the  accuracy  of  the  Gauss  method  dropped  off  sharply.  Gauss  suffered  from 
poor  estimates  of  the  slant  range  when  the  noise  level  in  the  data  was  as  large  as  .01 
degrees.  Over  data  sets  spanning  less  than  .5  degrees  of  arc  the  Gauss  technique  began 
to  fail  entirely  and  it  becomes  necessary  to  use  the  spherical  trigonometry  technique. 
The  spherical  trigonometry  technique  did  suffer  from  poor  results  when  looking  off 
nadir,  but  as  long  as  measurements  remained  within  10  degrees  of  nadir  the  least 
squares  process  was  easily  able  to  compensate. 

The  effects  of  ground  sample  distance  on  the  ability  of  these  methods  to  ef¬ 
fectively  determine  the  orbit  of  a  fastwalker  was  also  analyzed,  ft  was  found  that 
for  any  reasonable  value  the  level  of  noise  in  the  data  that  would  be  imparted  by  a 
large  ground  sample  distance  did  not  have  a  significant  affect  on  the  results.  Since 
commercially  available  platforms  operating  in  the  near  infrared  spectrum  are  avail¬ 
able  with  ground  sample  distances  of  10  meters  or  less  from  LEO,  it  is  reasonable  to 
assume  that  future  platforms  placed  in  a  HEO  or  GEO  orbit  will  not  have  ground 
sample  distances  much  larger  than  one  kilometer.  In  fact,  the  level  of  noise  imparted 
by  large  ground  sample  distance  does  not  prevent  identification  of  the  target  satellite 
until  these  values  are  greater  than  12.65  kilometers. 
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Non-linear  least  squares  proved  to  be  an  effective  method  of  refining  the  estimate 
of  the  state  obtained  from  the  initial  orbit  determination  process.  It  allowed  for  the 
use  of  additional  data  and  proved  robust  enough,  even  with  very  short  data  sets,  to 
correct  for  several  degrees  of  error  in  the  initial  guess  of  the  state.  Although  the 
Gauss  initial  orbit  determination  method  sometimes  provided  an  initial  guess  that 
was  very  accurate  over  data  sets  spanning  several  minutes,  its  accuracy  could  not  be 
confirmed  without  taking  into  account  the  additional  data  through  the  least  squares 
process.  This  made  the  non-linear  least  squares  process  necessary  to  confirm  the 
accuracy  of  the  state  vector  developed  using  Gauss.  When  the  target  was  in  a  highly 
eccentric  orbit,  the  method  did  not  fail  but  performed  poorly  even  when  the  target 
was  relatively  close  to  the  surface  of  the  Earth.  The  method  was  unsuccessful  in 
identifying  targets  at  or  near  GEO. 

Searching  for  a  fastwalker  target  based  on  its  right  ascension  of  the  ascending 
node  and  its  inclination  proved  to  be  an  effective  way  of  identifying  the  orbit.  This 
method  is  heavily  dependent  on  using  an  Air  Force  Space  Command  catalog  published 
during  the  period  that  the  fastwalker  was  observed.  While  inclination  is  nearly  con¬ 
stant,  the  node  is  not,  and  using  a  catalog  published  even  days  apart  could  invalidate 
the  method.  However,  these  catalogs  are  published  several  times  a  week  and  older 
copies  are  readily  available.  The  argument  of  latitude  was  effective  as  the  final  step 
in  confirming  the  identity  of  the  fastwalker  and  even  at  the  most  crowded  inclinations 
there  was  no  difficulty  in  identifying  the  target. 

The  method  developed  was  not  intended  to  determine  the  orbit  of  the  target 
alone.  It  was  effective  in  identifying  the  target  fastwalker  from  the  data  in  the  Space 
Command  catalog.  In  not  attempting  to  fully  determine  the  orbital  elements  of  the 
target  fastwalker  directly,  this  method  circumvented  some  of  the  problems  which  could 
have  arisen.  Over  very  short  data  tracks  with  only  angular  data,  an  accurate  estimate 
of  the  orbital  elements  of  the  target  would  have  been  impossible.  In  the  recent  paper 
by  the  Chinese  National  Defense  Institute  [7]  over  30  minutes  of  data  was  needed 
in  order  to  obtain  even  a  poor  independent  estimate  of  the  target  satellite  from  a 
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HEO  or  GEO  orbit,  excluding  the  fact  that  it  is  unlikely  that  a  LEO  satellite  would 
be  visible  to  a  satellite  in  HEO  orbit  for  that  amount  of  time.  Using  much  shorter 
data  spans,  it  was  possible  to  take  advantage  of  the  existence  and  accuracy  of  the 
Space  Command  catalog  and  using  only  a  few  orbital  elements  identify  the  correct 
satellite  in  the  catalog  when  the  orbit  is  determined  to  be  near  LEO.  From  the  two 
line  element  set  in  the  catalog  the  orbit  is  then  known  to  an  acceptable  standard. 

5. 1  Recommendations 

The  topic  of  fastwalker  orbit  determination  has  not  received  much  attention  in 
recent  years.  It  was  determined  to  be  impossible  for  single  data  tracks  in  a  geosyn¬ 
chronous  orbit  and  no  other  orbit  types  were  investigated.  The  advantage  in  working 
in  a  HEO  orbit  was  the  fact  that  these  orbits  are  not  crowded  by  other  platforms,  so 
the  key  assumption  that  the  object  was  located  near  the  Earth’s  surface  was  made. 
In  Osedacz’s  thesis  [9],  he  determined  that  the  majority  of  fastwalkers  visible  to  the 
platform  he  was  working  with  were  in  fact  also  in  geosynchronous  orbit.  This  led  to 
his  inability  to  determine  the  range  between  the  sensor  and  the  target. 

The  focus  of  this  thesis  was  on  determining  the  orbit  of  the  target  fastwalker. 
Very  little  attention  was  given  to  the  actual  nature  of  the  fastwalker  phenomenon. 
Fastwalkers  represent  objects  large  enough  to  reflect  enough  light  in  one  direction  to 
be  clearly  visible  from  a  HEO  or  GEO  orbit.  When  these  events  are  likely  to  occur 
and  under  what  circumstances  needs  further  analysis.  In  order  to  do  so,  significant 
amounts  of  data  on  actual  fastwalker  observations  would  be  necessary.  The  time  of 
the  observation  and  the  location  of  the  observer  would  be  the  needed  to  even  begin  a 
thorough  analysis.  Since  this  data  was  not  available  for  this  thesis  the  topic  was  only 
briefly  discussed.  This  is  an  area  where  a  significant  amount  of  future  work  could  be 
conducted. 

There  are  several  areas  where  future  work  could  be  conducted  to  further  the 
results  obtained  here.  Obtaining  actual  fastwalker  data  and  attempting  to  conduct 
a  successful  orbit  matching  using  the  techniques  developed  in  this  thesis  would  be  a 
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useful  follow  on  project.  Additionally,  future  work  could  be  conducted  by  attempting 
to  use  the  same  orbit  matching  idea  and  developing  a  matching  technique  which 
accounts  for  model  error.  This  would  improve  the  values  of  the  covariance  matrix. 
If  the  accuracy  of  these  values  is  improved  than  the  entire  process  from  data  entry 
into  the  system  until  target  identification  could  be  automated  removing  the  need  for 
a  user  search  of  the  catalog. 

During  the  development  of  this  thesis,  I  was  also  impressed  by  the  almost  com¬ 
plete  lack  of  modern  research  into  angles-only  initial  orbit  determination  techniques. 
Gooding,  [4, 5]  is  the  only  person  currently  active  in  researching  new  techniques,  al¬ 
though  Chris  Sabol  [11,12]  is  working  on  applications  in  improving  the  accuracy  of 
entries  in  the  Air  Force  Space  Command  catalog  with  angular  data.  For  my  purposes, 
I  found  the  Gauss  method  when  used  in  conjunction  with  the  Herrick-Gibbs  method 
provided  an  estimate  of  the  location  of  the  orbit  which  was  more  than  adequate.  Ad¬ 
ditional  experimentation  on  refining  the  least  squares  process  is  needed  to  attempt  to 
account  for  model  error  in  the  covariance  matrix.  With  model  error  accounted  for,  it 
would  be  possible  to  automate  the  entire  process  and  remove  the  manual  searching 
of  the  catalog  currently  necessary  by  the  user  to  characterize  the  target.  Despite  the 
success  of  this  method,  determining  the  exact  orbit  of  another  satellite  in  GEO  from 
a  GEO  sensor  using  angles-only  data  may  still  remain  impossible.  As  it  currently 
stands  this  method  allows  for  the  identification  of  over  82  percent  of  the  objects  listed 
in  the  AFSPC  catalog. 
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Appendix  A.  Residual  Plots 

This  Appendix  includes  all  of  the  residual  charts  referenced  in  chapter  4.  The  residual 
plots  shown  are  all  with  10e-7  degrees  of  error.  For  examples  of  residual  plots  with 
larger  errors  please  refer  to  chapter  4. 


Azimuth  Residuals 


Figure  A.l:  Azimuth  Residuals  With  3.18e-3  Noise  for  Cosmos 
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Elevation  Residuals 
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Figure  A. 2:  Elevation  Residuals  With  3.18e-3  Noise  for  Cosmos 


Figure  A. 3:  Azimuth  Residuals  With  3.18e-3  Noise  for  SL-8  R/B 
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Elevation  Residuals 


Figure  A. 4:  Elevation  Residuals  With  3.18e-3  Noise  for  SL-8  R/B 


Azimuth  Residuals 


Figure  A. 5:  Azimuth  Residuals  With  3.18e-3  Noise  for  Jason 
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Elevation  Residuals 


Figure  A. 6:  Elevation  Residuals  With  3.18e-3  Noise  for  Jason 


Azimuth  Residuals 


Observation  # 

Figure  A. 7:  Azimuth  Residuals  With  3.18e-3  Noise  for  1SS 
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Elevation  Residuals 


Figure  A. 8:  Elevation  Residuals  With  3.18e-3  Noise  for  ISS 
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Figure  A. 9:  Azimuth  Residuals  With  3.18e-3  Noise  for  N2  RB 
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Elevation  Residuals 


Figure  A. 10:  Elevation  Residuals  With  3.18e-3  Noise  for  N2  RB 


Figure  A. 11:  Azimuth  Residuals  With  3.18e-3  Noise  for  ISS 
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Elevation  Residuals 


Figure  A. 12:  Elevation  Residuals  With  3.18e-3  Noise  for  ISS 
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Figure  A. 13:  Azimuth  Residuals  With  3.18e-3  Noise  for  OPS  4467  A 
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Figure  A.  14:  Elevation  Residuals  With  3.18e-3  Noise  for  OPS  4467  A 
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Figure  A. 15:  Azimuth  Residuals  With  3.18e-3  Noise  for  Delta  2  R/B 
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Figure  A. 16:  Elevation  Residuals  With  3.18e-3  Noise  for  Delta  2  R/B 
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Appendix  B.  Matlab  core  code 

This  appendix  contains  the  core  code  elements.  All  code  was  written  and  functions 
in  Matlab®  . 


B.l  Gauss  Initial  Orbit  Determination 

Listing  B.l:  The  Gaussian  initial  orbit  determination  code. 
(appendix2  /  GaussIOD .  m) 

17.  ... 


7. 

% 

7. 

7.  this  code  is 

6  % 

7.  inputs 
7.  year 
7.  month 
7.  day 

11  7.  hour  1-3 
7.  si -3 
7.  alpha 
7,  sigma 
7.  r_sitel  -3 
16  7. 

7.  outputs 

7. 

% 

7. 

21  7.  ... 


Gauss_I0D . m 

the  Gauss  initial  orbit  determination 


description 
Of  observation 

II  II 


range  /  units 
Time 


3  hours  of  observation 
3  seconds  of  observation 
Azimuth 
Elevation 

location  of  sensor  ECI 


radians 

radians 

kms 


r2 , v2  , 

argument  of  lat  , omega  (RAAN) 
i  (inclination), 

Thetadot  (alt) 


km/  km/s 
degrees 

NLLS  input  term 


clc 

clear  all 
close  all 
26  /.Constants 

u_earthkm3per se  c2 =398600 . 4415 ; 
radius_earth=6378 . 1363 

7,  for  Cosmos 

31  “/.Time  of  Observations 

year=2007 ; month=9; day =27 ; 
hourl=15 ; hour2=15; hour3=15; 
mini =48 ; min2=49 ; min3=51 ; 
sl=0.74;s2=59.74;s3=0.74; 

36  °/0R_site 
7.  data 
7.  Azimuth 
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alpha = [0 . 601631211134485 ; 1.057947 399382605; 1.306829387830625]; 
“/oElevat  ion 


41  sigma =[-1.49 83077489 6 3798; -1.506215139106105; -1.504976881780819]; 
r_sitel=[-17470. 7644479361  12124.4096188952 

27495 . 1031194714] ’ ; 

r_site2= [-17582 . 718937935  11991.5070697076 

27779 . 4780612869]  ’  ; 

r_site3= [-17639 . 2080231445  11922.7654824776 

27923 . 8294989853] ’ ; 


46  '/.Line  of  Sight  Unit  Vector 

Ll  =  [-  cos ( sigma (1 ,1) )*cos (alpha (1 ,1) )  ; cos (sigma(l ,1) )*sin(alpha.  .  . 
(1  ,  1)  )  ; sin (sigma  (1,1))]; 

L2=[-cos( sigma (2 , 1) ) *cos (alpha (2 , 1) )  ; cos (sigma (2 , 1) ) *sin (alpha.  .  . 

(2.1) ) ; sin (sigma (2,1))] ; 

L3=[-cos( sigma (3,l))*cos(alpha(3,l));cos(sigma(3,l))*sin(alpha... 

(3.1) )  ; sin ( sigma  (3,1))]  ; 


51  r_site=[r_sitel ,r_site2 ,r_site3] ; 

7. - 


56 


61 


66 


7.  Quadrant  check  of 
r_sites  are 
°/0in  same  quadrant) 
if  r_sitel(l,l)  <  0 
display ( ’ Sensor 
modifier=l ; 
elseif  r_sitel(l,l) 
display ( ’ Sensor 
modif i er  =0 ; 
elseif  r_sitel(l,l) 
display ( ’ Sensor 
modif i er  =0 ; 
elseif  r_sitel(l,l) 
display ( ’ Sensor 
modifier=l ; 


sensor  platform  in  ECI  frame  (assumeing  all 


&&  r_sitel(2,l)  <  0 
in  quadrant  3 ’ ) 

>  0  &&  r_sitel(2,l)  >  0 
in  quadrant  1 ’ ) 

<  0  &&  r_sitel(2,l)  >  0 
in  quadrant  2 ’ ) 

>  0  &&  r_sitel(2,l)  <  0 
in  quadrant  4 ’ ) 


end 


7 - 

7.  CONVERSION  OF  AZ  EL  TO  RIGHT  ASCENSION  AND  DECLINATION 
71  if  modifier  ==  0 
x=[l ,0,0] ’ ; 

r _s it e lxy  = [ r _sit e 1  (  1  ,  1 )  ,r_sitel(2,l)  ,0]  ’; 

phi  =  acos (dot (r_sitelxy , r_sitel) /( (norm (r_sitelxy) *norm(r_sitel) ) ) ) .  .  . 

f 

76  theta=acos (dot (r_sitelxy , x) / ( (norm (r_sitelxy) *norm (x) ) ) ) ; 

SEZ= [sin (phi ) *cos (theta)  , -sin (theta)  , cos ( theta )* cos ( phi )  ;  .  .  . 
sin (phi ) *sin (theta)  , cos (theta)  , s in ( thet a ) *  cos ( phi )  ;  .  .  . 

-  cos ( phi )  ,0,sin(phi)]  ; 

L 1  =  SEZ  *  L 1  ; 

81  r_site2xy=[r_site2(l,l)  ,r_site2(2,l)  ,0]  ’; 
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phi=acos (dot (r_site2xy , r_site2) / ( (norm (r_site2xy) *  norm (r_site2) ) ) ) 

i 

theta  =  acos (dot (r_site2xy , x) / ( ( norm (r_site2xy) *  norm (x ) ) ) )  ; 

SEZ =[ sin ( phi ) *  cos ( theta)  , -sin (theta)  , cos (theta) *cos (phi)  ;  .  .  . 

86  sin ( phi ) * s in ( theta )  , cos (theta)  , s in ( thet a ) *  cos ( phi )  ;  .  .  . 

-  cos ( phi )  ,0,sin(phi)]  ; 

L2  =  SEZ  *  L2  ; 

r_site3xy=[r_site3(l,l)  ,r_site3(2,l)  ,0]  ’  ; 

phi=acos (dot (r_site3xy , r_site3) / ( (norm (r_site3xy) *  norm (r_site3) ) ) ) 
> 

91  theta=acos (dot (r_site3xy , x) / ( ( norm (r_site3xy) *norm (x ) ) ) ) ; 

SEZ =[ sin ( phi ) *  cos ( theta)  , -sin (theta)  , cos ( thet a ) *  cos ( phi )  ;  .  .  . 
sin (phi ) *sin (theta)  , cos (theta)  , s in ( thet a ) *  cos ( phi )  ;  .  .  . 

-  cos ( phi )  ,0,sin(phi)]  ; 

L3  =  SEZ  *  L3  ; 

96  elseif  modifier  ==  1 
x=[l  ,0 ,0]  ’  ; 

r_sitelxy=[r_sitel(l,l)  ,r_sitel(2,l)  ,0]  ’ ; 

phi=acos (dot (r_sitelxy , r_sitel) /( (norm (r_sitelxy) *norm(r_sitel) ) ) ) 

f 

101  theta=2*pi-acos(dot(r_sitelxy ,x)/((norm(r_sitelxy)*norm(x)))) ; 

SEZ = [ sin ( phi ) *  cos ( thet a)  , -sin (theta)  , cos ( thet a ) *  cos ( phi )  ;  .  .  . 
sin (phi ) *sin (theta)  , cos (theta)  , s in ( thet a ) *  cos ( phi )  ;  .  .  . 

-  cos ( phi )  ,0,sin(phi)]  ; 

L 1  =  SEZ  *  L 1  ; 

106  r_site2xy=[r_site2(l,l)  ,r_site2(2,l)  ,0]  ’ ; 

phi=acos (dot (r_site2xy , r_site2) / ( (norm (r_site2xy) *  norm (r_site2) ) ) ) 
) 

theta=2*pi-acos (dot (r_site2xy ,x) /( (norm(r_site2xy) *  norm (x) ) ) ) ; 

SEZ= [sin (phi ) *cos (theta)  , -sin (theta)  , cos ( thet a ) *  cos ( phi )  ;  .  .  . 

Ill  sin (phi ) *sin (theta)  , cos (theta)  , s in ( thet a ) *  cos ( phi )  ;  .  .  . 

-  cos ( phi )  ,0,sin(phi)]  ; 

L2  =  SEZ  *  L2  ; 

r_site3xy=[r_site3(l,l)  ,r_site3(2,l)  ,0] 

%x=[r_site3(l,l)  ,0,0]  ’  ; 

116  phi=acos (dot (r_site3xy , r_site3) / ( (norm (r_site3xy) *norm (r_site3) ) ) ) 

i 

theta=2*pi-acos (dot (r_site3xy ,x) /( (norm(r_site3xy) *norm (x) ) ) ) ; 

SEZ= [sin (phi ) *cos (theta)  , -sin (theta)  , cos ( thet a ) *  cos ( phi )  ;  .  .  . 
sin (phi ) *sin (theta) , cos (theta) , sin (theta) *cos (phi ) ; . . . 

-  cos ( phi )  ,0,sin(phi)]  ; 

121  L3  =  SEZ  *  L3  ; 

end 

% - 

L  =  [LI  ,  L2  ,  L3]  ; 

“/-Figure  out  angles - 
126  dec=asin (L3 (3 , 1) ) *180/pi ; 
csigmat  =  sqrt (1-L3 (3 , 1) )  ; 
ra=asin(L3 (2 , 1) / csigmat ) *180/pi ; 
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”/o  Beginning  Gaussian  Method 

131  t aul = ( hour 1 *60+minl + s 1 /60) - (hour2 *60+min2+s2 /60) ; 
t au3 = ( hour3 *60+min3 + s3 /60) - (hour2 *60+min2+s2 /60) ; 
taul_sec=taul*60; 
tau3_sec=tau3*60; 

136  al=tau3/ (tau3 -taul ) ; 
a3=-taul / (tau3 -taul ) ; 

•/. - 

7,  convert  time  to  seconds 

alu=(tau3_sec*((tau3_sec-taul_sec) ~2-tau3_sec~2) ) / (6*(tau3_sec-. .  . 
taul.sec) ) ; 

141  a3u=-(taul_sec*((tau3_sec-taul_sec)~2-taul_sec~2))/(6*(tau3_sec-.. 
taul.sec) ) ; 

invL=inv (L) ; 

M=invL*r_site ; 

dl=M(2,l)*al-M(2,2)+M(2,3)*a3; 

146  d2=M (2 , 1) *alu+M(2 ,3) *  a3u  ; 

C=L2 ’*(r_site2) ; 

70Calculate  the  8th  order  polynomial’s  roots 

poly=  [1  0  -(dl ~2+2*C*dl+(r_site2  ’*r_site2) )  0  0  (-2*... 

u_earthkm3per se c2 * ( C* d2+dl * d2 ) )  0  0  -u_earthkm3persec2 ~2*d2 ~2] ; 
r2=roots (poly ) ; 

151  format  long 
display (r2) 

r2= input ( 1  Enter  the  correct  root  choice  numerically:  ’); 

ul  =  u_e  art hkm3per s  e  c2 / r2 " 3 ; 
c=[al+alu*ul ; -1; a3+a3u*ul] ; 

156  p= inv ( [ c ( 1  ,  1 )  0  0;0  c(2,l)  0;0  0  c (3 , 1) ] ) *M*-c ; 

pcheck  =  inv ( [c (1  ,  1)  0  0 ; 0  c(2,l)  0;0  0  c (3 , 1) ] ) *M*-c ; 

magp  =  ( sqrt (p(l,l)"2  +  p(2,l)~2  +  p(3,l)~2))  ; 

7»we  now  have  first  guess  at  position  vectors  rl,r2,r3 
161  r=[p(l,l)*Ll+r_sitel ,p(2,l)*L2+r_site2,p(3,l)*L3+r_site3] ; 
rl=r  (  :  ,1)  ; 
r2  =  r  (  :  ,2)  ; 
r3  =  r  (  :  ,3)  ; 

166 

70We  now  use  either  the  Gibbs  or  Herr  i  ck  -  Gibbs  method  to  calculate 
the  velocity  vector 

alphal  =  acos ((rl  ’  *  r2  ) / ((sqrt (r2 (1 , 1) ~  2  +  r2 (2 , 1) ~2+r2  (3 , 1) ~ 2 )  )  * ( sqrt  ( 
rl (1 , 1) "2+rl (2 , 1) "2+rl  (3 , 1) ~2) ) ) ) ; 
alpha3=acos ((r3 ’  *r2) / ((sqrt (r2 (1 , 1) ~2+r2 (2 , 1) ~2+r2 (3 , 1) ~2 ) ) *(sqrt ( 
r3(l  ,1) "2+r3 (2 , 1) ~2+r3(3 ,1) ~2) ))) ; 
alpha_total=alphal *180/pi+alpha3*180/pi ; 

171 

if  alpha_total  <  5  7»Herrick -Gibbs 

tau31=(hour3*3600+min3*60+s3) -(hourl*3600+minl*60+sl) ; 

tau21=(hour2*3600+min2*60+s2) -(hourl*3600+minl*60+sl) ; 

Z23  =  cross(r2,r3)  ; 
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176 


magr 1 = ( sqrt (rl(l,l)~2+rl(2,l)~2+rl(3,l)~2)); 
magr 2  =  sqrt (r2 (1 , 1) ‘2  +  r2 (2 , 1) ~2  +  r2  (3 , 1) ~2)  ; 
magr3  =  (sqrt (r3(l ,1) ~  2  +  r3 (2 ,1) ~2+r3(3,l) ~  2 ) )  ; 

'/.coplanar  check 

alphacoplanar=90*pi/180-acos ((Z23 ’ * r 1 ) / ((sqrt (Z23(l  ,1) ~  2  +  Z23 . . 
(2,1) ~  2  +  Z23 (3 ,1) ~2) ) *(sqrt (rl(l,l)~2+rl(2,l)~2+rl(3,l)~2))) 
)  ; 

181  if  alphacoplanar  >  .001 

display (’ Vectors  are  not  Coplanar! ’) 
break 

end 

v2=-tau3*(l/ (tau21*tau31) +u_e  art hkm3per  s  ec2 /(12*magrl~3))*rl+( 
t au3 -tau21) *(l/(tau21*tau3) +u_ear thkm3per se  c2 / (12*magr2~3) ) 
*r2+tau21*(l/(tau3*tau31) +u_earthkm3persec2 /(12*magr3~3) ) * . 
r3 

186 

elseif  alpha_total  >=  5  % Gibbs 
Z12=cross(rl ,r2) ; 

Z23  =  cross (r2 , r3)  ; 

Z31  =  cross(r3,rl)  ; 

191  '/.coplanar  check 

alphacoplanar=90*pi/180-acos ((Z23 ’ * r 1 ) / ((sqrt (Z23(l  ,1) ~  2  +  Z23 . . 
(2,1) ~  2  +  Z23 (3 , 1) ~2) ) * ( sqrt (rl(l,l)~2  +  rl(2,l)~2  +  rl(3,l)~2))) 
)  ; 

if  alphacoplanar  >  .001 

display (’ Vectors  are  not  Coplanar!  ’) 
break 

196  end 

N=( (sqrt (rl(l,l)~2+rl(2,l)‘2+rl(3,l)~2) ) *Z23  +  (sqrt (r2 (1  ,  1) ~2+ . 
r2 (2 , 1) "2+r2 (3 ,1) ~2) ) *Z31+(sqrt (r3 (1 , 1) “2+r3 (2 , 1) ~2+r3 (3 , 1) 
~2)  )  *Z12)  ;  7.ER-3 

magr 1 = ( sqrt (rl(l,l)‘2+rl(2,l)~2+rl(3,l)“2)); 
magr 2  =  sqrt (r2 (1 , 1) ‘2  +  r2 (2 , 1) ~2  +  r2  (3 , 1) ~2)  ; 
magr3  =  (sqrt (r3(l ,1) ~2  +  r3(2 ,1) ~2+r3(3,l) ~2 ) )  ; 

201  D=Z12+Z23+Z31 ; 

S= ( magr 2 -magr3)*rl+(magr3-magrl)*r2+( magr 1 -magr2) *r3 ; 

B=  cross (D , r2 )  ; 

Lg  =  sqrt (u_earthkm3persec2 / (N  ’  *D) )  ; 
v2  =  (  Lg/magr2  )  *B  +  Lg  *  S  ;  0/0ER/TU 

206  end 

“/.needed  orbital  elements  from  r  and  v  to  get  pi 

h=cross  (r2  ,  v2)  ;  “/angular  momentum 

magv2  =  ( sqrt (v2 (1 , 1) ~2  +  v2 (2 , 1) ~2  +  v2  (3 , 1) ~2) ) ; 

211  s igma= ((magv2~2)/2)-( u_earthkm3per se  c2 /magr 2 )  ; 

e  =  1 / u_earthkm3per sec2  * ( (magv2  ~2 - ( u_earthkm3per se  c2 /magr2) ) *r2-(r2 . 

’ *v2)  1  *v2)  ; 

mage  =  (sqrt (e(l  ,1) ~2+e(2,l) ~2+e(3,l) ~2) )  ; 
a  =  l ; 
pl  =  1 ; 

216  if  mage  <  1 

a=-u_earthkm3persec2 / (2*sigma) ; 
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pl=a* (1 -mage ~2) ; 
elseif  e  >=  1 

display (  ’  Orbit  is  parabolic  with  eccentricity  equal  to  or 
greater  than  1’); 

221  end 


'/.now  that  we  can  determine  the  semi -parameter  pi  we  can  get  a  new 
estimate 

226  '/,of  the  slant  range 

deltavl=-acos ((rl/magrl) ,*r2/magr2) ; 
deltav3=acos ((r3/magr3) ,*r2/magr2) ; 
vll=acos((e’*rl)/(mage*magrl)) ; 
v22=acos((e’*r2)/(mage*magr2)) ; 

231  v33=acos((e’*r3)/(mage*magr3)) ; 

fl=l-(magrl/pl) *(l-cos (deltavl) ) ; 
f3=l-(magr3/pl) * ( 1 -  cos (deltav3) ) ; 

gl  =  ( (magrl ) * (magr2) *  s  in (deltavl ) ) / sqrt ( u_ear thkm3per se  c2  *pl )  ; 
g3  =  ( (magr3) * (magr2) *  s  in (deltav3) ) / sqrt ( u_ear thkm3per se  c2  *pl )  ; 

236  cl=g3/(fl*g3-f3*gl) ; 
c3=-gl/(fl*g3-f3*gl) ; 

p=inv([cl  0  0;0  c(2,l)  0;0  0  c3])*M*-c; 

pcheckl = inv ( [ c 1  0  0;0  c(2,l)  0;0  0  c3])*M*-c 

magpl  =  ( sqrt (p ( 1  ,1) ~2+p(2,l) ~2  +  p(3,l) ~2 ) )  ; 

241 


'/.now  that  we  have  two  guesses  at  slant  range  we  need  to  iterate  to 
ref ine 

'/the  initial  estimate  of  slant  ranges 
count  =  1 ; 

246  while  abs (pcheck (1 , 1) -pcheckl  (1  ,  1)  )  >  .1  &&  abs ( pcheck (3 , 1 ) - . .  . 
pcheckl (3,1))  >  .1 

pcheckl =pcheck ; 

r=[p(l ,l)*Ll+r_sitel ,p(2,l)*L2+r_site2 ,p(3,l)*L3+r_site3] ; 
r 1 =r  (  :  ,1)  ; 
r2  =  r  (  :  ,2)  ; 

251  r3=r(: ,3) ; 

°/«We  now  use  either  the  Gibbs  or  Herr  i  ck  -  Gibbs  method  to  ... 
calculate  the  velocity  vector 

alphal  =  acos((rl ’*r2)/((sqrt(r2(l  ,1) ~2+r2(2,l) ~2+r2(3,l) ~2) ) *  (  . 

sqrt (rl (1 , 1) “2+rl (2 , 1) ~2+rl  (3 , 1) ~2) ) ) )  ; 
alpha3  =  acos ( (r3 ’  *r2) /( (sqrt (r2 (1 , 1) ~2+r2 (2 , 1) ~2  +  r2 (3 , 1) ~2) ) * ( . 

sqrt (r3 (1 , 1) ~2+r3 (2 , 1) "2+r3 (3 , 1) ~2) ) ) ) ; 
alpha_total=alphal *180/ pi+alpha3*180/pi ; 

256 

if  alpha_total  <  5  %Herri ck - Gibbs 

tau31=(hour3*3600+min3*60+s3) -(hourl*3600+minl*60+sl) ; 
tau21=(hour2*3600+min2*60+s2) -(hourl*3600+minl*60+sl) ; 
Z23=cross (r2 , r3) ; 

261  '/.coplanar  check 
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alphacoplanar=90*pi/180-acos((Z23’*rl)/((sqrt(Z23(l,l)~2+. 
Z23  (2,1) ~  2  +  Z23 (3 ,1) ~2) ) *(sqrt (rl(l,l)~2+rl(2,l)~2+rl... 

(3.1)  “2))))  ; 

if  alphacoplanar  >  .0001 

display (’ Vectors  are  not  Coplanar ! ' ) 
break 

266  end 

v2=-tau3*(l/ ( t au2 1  *  t au3 1 ) +u_eart hkm3per se  c2 / (12*magrl ~3) ) * 
rl  +  (tau3-tau21) *(1/ (tau21*tau3) +u_earthkm3per se  c2 / (12*. 
magr2~3))*r2+tau21*(l/(tau3*tau31) +u_earthkm3persec2 . .  . 
/(12*magr3~3) )*r3 

elseif  alpha_total  >=  5  "/.Gibbs 
Z12=cross (rl , r2) ; 

271  Z23=cross (r2 , r3) ; 

Z31=cross (r3 , rl) ; 

"/.coplanar  check 

alphacoplanar=90*pi/180-acos((Z23’*rl)/((sqrt(Z23(l,l)~2+. 
Z23 (2,1) ~  2  +  Z23 (3 , 1) ~2) ) *( sqrt (rl(l,l)~2+rl(2,l)~2+rl... 

(3.1) -2))))  ; 

if  alphacoplanar  >  .0001 

276  display (’ Vectors  are  not  Coplanar!") 

break 

end 

N=( (sqrt (rl (1 ,1) ~2+rl (2 ,1) ~2  +  rl (3,1) ~2) ) *Z23+(sqrt (r2 (1 , 1) 
~2+r2(2,l) ~2+r2(3,l) ~2) ) *Z31+(sqrt (r3 (1 , 1) ~2+r3 (2 , 1) ~2+ 
r3  (3 , 1)  ‘2)  )  *Z12)  ;  "/,ER~3 
magr 1 = ( sqrt (rl(l,l)~2+rl(2,l)~2+rl(3,l)~2)); 

281  magr 2  =  sqrt (r2 (1 , 1) "2+r2 (2 , 1) ~2+r2 (3 , 1) ~2)  ; 

magr3  =  (sqrt (r3 (1 , 1) ~2+r3 (2 , 1) ~2  +  r3 (3 , 1) ~2) )  ; 

D=Z12+Z23+Z31 ; 

S=(magr2-magr3)*rl+(magr3-magrl)*r2+(magrl-magr2)*r3; 

B  =  cross (D , r2)  ; 

286  Lg=sqrt (u_earthkm3persec2/ (N  ’  *D) )  ; 

v2  =  (  Lg/magr  2  )  *B  +  Lg*S  ;  7.ER/TU 

end 

"/.needed  orbital  elements  from  r  and  v  to  get  pi 
h=  cross  ( r2  ,  v2  )  ;  "/.angular  momentum 
291  magv2  =  (  sqrt  (v2  (1  ,  1)  '2+v2  (2 , 1)  ~2  +  v2  (3 , 1)  ~2)  )  ; 

sigma=((magv2~2) / 2) - ( u_earthkm3per se  c2 /magr 2 )  ; 

e  =  1 / u_eart hkm3per se c2  * ( ( magv2  ~2 - ( u_earthkm3per sec2 /magr2) ) *r2 . 
- ( r2 ’ * v2 )  ’*v2)  ; 

mage  =  (sqrt(e(l ,1) ~2+e(2,l) ~2  +  e(3,l) ~2) )  ; 
a  =  l ; 

296  pi  =  1 ; 

if  mage  <  1 

a=-u_earthkm3persec2 / (2*sigma) ; 
pi =a* ( 1 -mage ~2) ; 
elseif  e  >=  1 

301  display (’ Orbit  is  parabolic  with  eccentricity  equal  to  or 

greater  than  1’); 

end 
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'/.now  that  we  can  determine  the  semi -parameter  pi  we  can  get  a 
new  estimate 
*/,of  the  slant  range 

delt avl =- acos ( ( r 1 /magr 1 ) ’*r2/magr2) ; 

306  deltav3=acos((r3/magr3) ’*r2/magr2) ; 

vll=acos((e’*rl)/(mage*magrl)) ; 
v22=acos((e’*r2)/(mage*magr2)) ; 
v33=acos((e’*r3)/(mage*magr3)) ; 
fl=l-(magrl/pl) * ( 1  -  cos (delt avl) ) ; 

311  f3=l-(magr3/pl)*(l-cos(deltav3)) ; 

gl=((magrl)*(magr2)*sin(deltavl))/sqrt( u_e  arthkm3per sec2  *pl )  ; 
g3  =  ( (magr3) * (magr2) *  sin (deltav3) ) /sqrt ( u_e  arthkm3per s  ec2  *pl )  ; 
cl=g3/(fl*g3-f3*gl) ; 
c3=-gl/(fl*g3-f3*gl) ; 

316  p=inv([cl  0  0;0  c(2,l)  0;0  0  c3])*M*-c; 

pcheckl=inv ( [cl  0  0 ; 0  c(2,l)  0;0  0  c3])*M*-c 

magpl  =  ( sqrt (p ( 1 ,1) ~2  +  p  (2  ,1) ~2  +  p (3 , 1 ) ~2)  )  ; 
count=count+l ; 

end 

321  [p  ,  a  ,  ecc  ,  incl  ,  omega  ,  argp  ,  nu  ,  m  ,  m_arglat  ]  =  rv2coe  (r2,v2); 

incl=incl*180/pi 
omega  =  omega *180  / pi 
m=m*180/pi 
argp  =  argp  *  1 80/  pi 

326  m_arglat =m_arglat *  180/pi 
'/.calc  estimated  alt 
h_target=norm(r2) 

thetadot  =  sqrt ( u_earthkm3per se  c2 /(h_target) ~  3) 

B.2  Spherical  Trigonometry  Initial  Orbit  Determination 

Listing  B.2:  The  spherical  trigonometry  initial  orbit  determination  code. 
(appendix2/SphericaltrigIODm.m) 

17.  ... 


7. 

7,  Spherical_trig_IOD  .  m 

7. 

°/0  this  code  is  the  initial  orbit  determination  method  using  lat 
long  data . 

6  % 

7.  inputs  description 

7,  day  =  dataimport  (  :  ,  1 )  ; 

7.  month  =  dat  aimport  (  :  ,  2)  ; 

7.  year  =  dat  aimport  (  :  ,3)  ; 

11  7.  hour  =  dataimport  (  :  ,4)  ; 

“/,  min  =  dataimport  (  :  ,5)  ; 

%  sec=dataimport(:,6); 

7.  Lat  =  d  at  aimport  (  :  ,  7)  ; 

%  Long  =  dataimport ( :  ,8)  ; 

16  7.  CASE1  :  Circular  ascending 


range  /  units 

TIME 

GEODEDIC 

SYMBOL 
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7. 

7. 

7. 

7. 

21  7. 
7. 
7. 
% 
7. 


CASE2 

CASE3 

CASE4 

outputs 


26  7. need  JDAY  .  M 

7.  ••• 


Circular  descending 
Polar  ascending 
Polar  desending 
: Omega  =  node 

1=  inclination 
Theta=Argument  of  latitude 
Thet adot = change  in  theta 
INITIAL  STATE 
[Omega  I  Theta  Thetadot] 

GSTIME.M,  MODULO. M 


clc 

format  long 

31 

fid  =  fopen(’ Thor RB_latlong.txt’  ,  ’  r  ’  )  ; 

dataimport  =  fscanf(fid,  ’7.g  7.g  7.g  7.g  7.g  7.g  7.g  7.g  ’  ,  [8  inf]) 
dat  aimport  =  dat  aimport  ’  ;  7.  It  has  two  rows  now. 
fclose (fid) 

36  day=dataimport ( : ,1) ; 
month  =  d at aimport (:  ,2)  ; 
year  =  dataimport ( :  ,3)  ; 
hour  =  dataimport ( :  ,4)  ; 
min  =  dat aimport ( :  ,5)  ; 

41  sec  =  dat aimport ( :  ,6)  ; 

Lat  =  dataimport(:  ,7)  ; 

Long  =  dataimport ( :  ,8)  ; 

46  for  i=l : length (Lat ) 

degree_lat  =  Lat (i  ,  1)  ; 
degree_lat=degree_lat*pi/180; 
fl at ening =1/298. 257223563; 
u=degree_lat ; 

51  geocentric_lat= at an ((1-flat ening) ~2*tan(u)) ; 

geocentric_latl=geocentric_lat  *180/ pi ; 

Lat (i , 1) =geocentric_latl ; 

end 

56  7oThere  are  4  cases  of  data  that  we  can  consider  using  here, 
have 

7»CASE1  :  Circular  ascending 
7.CASE2  :  Circular  descending 
7.CASE3:  Polar 

61  7«  These  cases  are  defined  in  terms  of  the  track  crossing  the 
equator  in  the 
7.  northern  hemisphere 

i  =  l; 


We 
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66 

71 

76 

81 

86 

91 

96 

101 

106 

111 


A=zeros((length(Long)-l)  ,4)  ; 
b=  []  ; 

7. - CHOOSE  ORBIT  DATA  TYPE 

70syms  easel  case2  case3 
7,  ie  choice  =  casel; 
choice=3 

7. - 


if  clioice==l 

for  i=l :( length (Long) -1) 

deltalat=Lat (i+1 , 1) -Lat (i , 1) ; 

de It  along  =(Long(i+l,l)-Long(i,l))*cos(Lat(i,l)*pi/180); 

alpha=atan2 (deltalat ,deltalong) ; 

b (i , 5) =alpha ; 

az =90 - ( alpha  *  1 80/ pi )  ; 

b  (  i  ,  3)  =  az  ; 

az  =  az  *pi / 180 ; 

7  step  1  7«  Omega  calc  here 

delta=atan(sin(Lat (i , 1) *pi/180) *tan (az) ) ; 

JD  =  JDAY ( year ( i  ,  1 )  ,  month ( i , 1 ) ,  day(i,l),  hour(i,l),  min(i,l) 
,  sec ( i  ,  1 )  )  ; 

[gst]  =  GSTIME  ( JD)  ; 
b ( i , 1 ) =gst  *  180/pi ; 
b(i , 2)=delta*180/pi ; 

"/.need  JDAY.M,  GSTIME.  M,  MODULO. M 
Omega=Long (i , 1) +gst*180/ pi-delta*180/pi ; 

A ( i , 1 ) =  Omega ; 
b ( i , 4) =  Omega ; 

'/.step  2 

theta  =  acos (cos (delta) *cos (Lat (i , 1) *pi/180) )  ; 
theta=theta*180/ pi ; 

A ( i ,3) =theta ; 

7  s  t  e  p  3 

inclination=asin(sin(Lat (i , 1) *pi/180) / sin(theta*pi/180) ) ; 
inclination=inclination*180/pi ; 

A(i ,2)=inclination; 

7  s  t  e  p  4 

thetadot=sqrt ((deltalong*pi/180) ~  2* ( cos (Lat (i , 1) *pi/180) ) ~ 2+ ( . 
deltalat*pi/180)~2)/((hour(i+l,l)*3600+min(i+l,l)*60+sec(i. 
+1 ,1) ) -(hour(i ,1) *3600+min(i ,1) *60+sec(i ,1) )) ; 

A(i ,4)=thetadot ; 
i=i+l ; 

end 

elseif  choice==2 

for  i=l :( length (Long) -1) 
deltalat=Lat(i  +  l,l)-Lat(i,l)  ; 

deltalong=(Long (i  +  1 , 1) -Long (i , 1) ) *cos (Lat (i , 1) *pi  /  180)  ; 
alpha=atan2 (deltalat ,deltalong) ; 
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116 

121 

126 

131 

136 

141 

146 

151 

156 

161 


b(i , 5)=alpha*180/pi ; 
az =90 - ( alpha  *  1 80/ pi )  ; 
bet a  =  90+ ( alpha  *  180/ pi )  ; 
b ( i ,6) =beta ; 
b ( i , 3) =az ; 
az  =  az  *pi / 180 ; 

'/.step  1  '/.Omega  calc  here 

delta=atan(sin(Lat (i , 1) *pi/180) *tan(beta*pi/180) )  ; 
delta=pi -delta; 

JD  =  JDAY ( year ( i  ,  1 )  ,  month(i,l),  day(i,l),  hour(i,l),  min(i,l) 
,  sec ( i  ,  1 )  )  ; 

[gst]  =  GSTIME  ( JD)  ; 
b(i , 1) =gst  *  180/ pi ; 
b(i ,2)=delta*180/pi ; 

'/.need  JDAY.M,  GSTIME. M,  MODULO. M 
Omega=Long (i , 1) +gst*180/ pi-delta*180/pi ; 

A ( i , 1 ) =  Omega ; 
b ( i , 4) =  Omega ; 

'/.step  2 

theta_tilde  =  acos (cos (delta) *cos (Lat (i , 1) *pi/180) )  ; 
theta=theta_tilde*180/pi ; 

A ( i ,3) =theta ; 

%  s  t  e  p  3 

inclination=asin(sin(Lat  (i  ,  1)  *pi/180)  /  sin(theta_tilde)  )  ; 
inclination=inclination*180/pi ; 

A(i ,2)=inclination; 

'/.step  4 

thetadot=sqrt ((deltalong*pi/180) ~  2* ( cos (Lat (i , 1) *pi/180) ) ~  2  +  ( . 
deltalat*pi/180) ~2)/((hour(i+l  ,1) *3600+min(i+l ,l)*60+sec(i. 
+l,l))-(hour(i,l)*3600+min(i,l)*60+sec(i,l))); 

A(i ,4)=thetadot ; 

i=i+l ; 

end 

elseif  choice==  3 

for  i=l :( length (Long) -1) 

deltalat=Lat(i  +  l,l)-Lat(i,l)  ; 

deltalong=(Long (i  +  1 , 1) -Long (i , 1) ) *  cos (Lat (i  ,  1) *pi  /  180)  ; 

alpha=atan2 (deltalat ,deltalong) ; 

b ( i , 5) =  alpha ; 

az =90 - ( alpha  *  1 80/ pi )  ; 

b ( i , 3) =az ; 

az  =  az  *pi / 180 ; 

‘/.step  1  %  Omega  calc  here 

delta=atan(sin(Lat (i , 1) *pi/180) *tan(az) ) ; 

JD  =  JDAY ( year ( i , 1 ) ,  month ( i , 1 ) ,  day(i,l),  hour(i,l),  min(i,l) 
,  sec ( i  ,  1 )  )  ; 

[gst]  =  GSTIME  ( JD)  ; 
b ( i , 1 ) =gst  *  180/pi ; 
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b(i ,2)=delta*180/pi ; 

'/.need  JDAY.M,  GSTIME.M,  MODULO. M 

Omega = Long (i , 1) +gst*180/pi-delta*180/pi ; 

A ( i , 1 ) =  Omega ; 

166  b ( i , 4) = Omega ; 

'/.step  2 

theta=acos (cos (delta) *  cos (Lat (i , 1) *pi  /  180) )  ; 
theta=theta*180/ pi ; 

A ( i ,3) =theta ; 

171 

‘/.step  3 

inclination=asin(sin(Lat (i , 1) *pi/180) / sin(theta*pi/ 180) ) ; 
inclination=180- inclination*180/pi ; 

A(i ,2)=inclination; 

176 

%  s  t  e  p  4 

thetadot  =  sqrt ((deltalong*pi/180) ~  2* ( cos (Lat (i , 1) *pi/180) ) ~  2  +  ( .  .  . 
deltalat*pi/180)  ~2)/((hour(i+l ,1) *3600  +  min(i  +  l ,l)*60+sec(i. .  . 
+l,l))-(hour(i,l)*3600+min(i,l)*60+sec(i,l))); 

A(i ,4)=thetadot ; 
i=i+l ; 

181  end 

end'/,  end  if  loop  for  choice 
display(A(i-l,  :)) 

display ([ 1  The  epoch  time  is  ’  , num2str (hour (i-1 , 1) )  ,  ’  :  ’  , num2str ( .  .  . 

min(i-l,l)),  ’  :  ’  , num2st r (sec(i-l,l))  ,]) 


B.3  Non-Linear  Least  Squares 

Listing  B.3:  The  least  squares  code.(appendix2/NLleastsquares.m) 

7.  ... 


7. 

7,  Non_linear_leastsquares  .  m 

4  l 

7«  this  code  is  the  master  Non-Linear  Least  Squares 

7. 

7.  inputs  description  range  /  units 

7.  day  =  dataimport  (  :  ,  1 )  ; 

9  7«  month  =  dataimport  (  :  ,2)  ; 

7«  year  =  data  import  (  :  ,3)  ; 

7«  hour  =  dataimport  (  :  ,  4)  ;  TIME 

7.  min  =  dataimport  (  :  ,5)  ; 

7.  sec  =  dataimport(:,6); 

14  7.  Xo  =  INITIAL  STATE=  [  Omega  I  Theta  Thetadot] 

7.  EPOCH  TIME 

7«  outputs  :  Omega  =  node 

/  I=inclination 

7«  Thet  a  =  Argument  of  latitude 

19  7«  Thet adot  =  change  in  theta 
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•/,  INITIAL  STATE 

°/0  [Omega  I  Theta  Thetadot] 

1  FILES  REQUIRED  TO 

7„  RUN  ;  testdatasensor4  .  txt  ;  testdatatarget4  .  txt  ;  G  .  M  ;  r.target  .  M;  .  .  . 
cooridnate.tra 
24  1  nsf  orm  .  m 
1  ... 


close  all 
clear  all 
clc 

29  meanr =10 ; 

'/.Cosmos 

Xo  =  [358 . 9086913581474*pi / 180 , 82 . 491055652040558*pi .  .  . 

/180  ,121.9992275169298*pi/180  ,  .0009138421573219655]  ; 
hour_epoch=15; 
min_epoch=49 ; 

34  sec_epoch=59 . 74 ; 

Xo  =  Xo  ’  ; 

'/.Gather  Data 
39  "/o  Target 

fid  =  f open (’ Cosmos_p . txt ’  ,  ’r’); 

dataimport  =  fscanf(fid,  ’  g  %g  °/,g  %g  %g  °/„g  °/,g  °/„g  lg’  ,  [9  inf]); 
dat  aimport  =  dat  aimport  ’  ;  %  It  has  two  rows  now. 
fclose (fid) 

44  day=dataimport ( : ,1) ; 
month  =  d at aimport (:  ,2)  ; 
year  =  dataimport ( :  ,3)  ; 
hour  =  dataimport ( :  ,4)  ; 
min  =  dat aimport ( :  ,5)  ; 

49  sec  =  dat aimport  (  :  ,6)  ; 
rx  =  dat  aimport  (  :  ,7)  ; 
ry  =  dataimport(:  ,8)  ; 
rz  =  dataimport(:  ,9)  ; 

54  '/.Sensor 

fid  =  f open (’ Cosmos_sensor_p . txt  ’  ,  ’r’); 

dataimport  =  fscanf(fid,  ’  °/  g  %g  °/,g  %g  °/0g  °/.g  °/.g  °/.g  lg’  ,  [9  inf]); 
dat  aimport  =  dat  aimport  ’;  %  It  has  two  rows  now. 
fclose (fid) 

59  day=dataimport ( : ,1) ; 
month  =  d at aimport (:  ,2)  ; 
year  =  dataimport ( :  ,3)  ; 
hour  =  dataimport ( :  ,4)  ; 
min  =  dat aimport ( :  ,5)  ; 

64  sec=d at aimport (: ,6) ; 

Rx=dataimport(: ,7) ; 

Ry  =  dat aimport  (  :  ,8)  ; 

Rz  =  dat aimport  (  :  ,9)  ; 
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69  [Az_ob  ,  El_ob]  =G  (rx  ,  ry  ,  rz  ,  Rx  ,  Ry  ,  Rz)  ; 
max_count=200; 
count  =  1  ; 

res_norm_old= 100000 ; 
res_norm_hist  =  []  ; 

74  temp=size ( Az_ob) ; 
sigma  =  0.0000001; 
nvecl=sigma*randn (temp) ; 
nvec2=sigma*randn (temp) ; 

Az_ob=Az_ob+nvecl ; 

79  El_ob=El_ob+nvec2 ; 

Q  =  diag (sigma“2*ones (length (Az_ob) *2 , 1) )  ; 

Qinv=inv (Q) ; 

ANS  =2 

while  ANS  >  1 
84  count  =  count  + 1 

'/.Begin  Non  Linear  Least  Squares 
for  i=l : length (Az_ob) 

%  Compute  current  state  from  inital  state  &  delta_t  (time  ... 
since  epoch) 

89  delta_t  =  -(hour_epoch*3600  +  min_epoch*60+sec_epoch)  +  (hour(i ,1)  .  . 

*3600  +  min (i , 1) *60+sec (i  ,  1) )  ; 
phi = [ 1  0  0  0 ; 0  1  0  0;0  0  1  delta_t;0  0  0  1]; 

Xn  =  phi *  Xo ; 


‘/.Calculate  pseudo_phi  using  Xn 
94  Qmega  =  Xn (1  ,  1)  ; 

I  =  Xn  (2 ,1)  ; 

Theta  =  Xn (3,1)  ; 

Thetadot =Xn (4 , 1 )  ; 
t=delta_t ; 

99  mue=  398601;  °/0  Earth  Gravitational  Parameter  (km“3/sec“2) 

xl  =  Omega ; 
x2  =  I  ; 
x3=Theta ; 
x4=Thetadot ; 

104  pseudo_phi  =  [  (  mue /x4  “  2)  “  (  1  /3)  *  (  -  sin  (  xl  )  *  cos  (  x3  )  -  cos  (  xl  )  *  cos  (  . 

x2) *sin (x3) )  ,  .  .  . 

(mue /x4 “ 2) “(1/3) *sin(xl)*sin(x2)*sin(x3) , ( mue /x4  “  2 ) “(1/3) . 
*  .  .  . 

(-cos(xl)*sin(x3)-sin(xl)*cos(x2)*cos(x3)) ,  -2/3/( mue / x4 . . 

“2) “ (2/3)  .  .  . 

*(cos (xl) *  cos (x3) -sin (xl) *cos (x2) *sin(x3) ) *mue / x4  ~  3 ;  .  .  . 

(mue /x4“2) “ (1/3) *(cos (xl) *  cos (x3) -sin (xl) *cos ( . 
x2  )  *  sin  (  x3  )  )  ,  .  .  . 

109  -  (mue /x4  “ 2)  “  (  1 /3 )  *  cos  ( xl  )  *  sin  ( x2  )  *  s in  ( x3  )  ,  (mue 

/x4 “ 2 )  “  (1/3)  .  .  . 

*(-sin(xl)*sin(x3)+cos(xl)*cos(x2)*cos(x3))  ,  .  .  . 
-2/3/ ( mue /x4“2) “(2/3) *(sin(xl)*cos(x3)+cos(xl) . 
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114 


119 


124 


129 


134 


139 


144 


*cos(x2)*sin(x3)) *mue/x4  ~3 ;  .  .  . 

0 ,  ( mue /x4  ~  2 )  ~(l/3)*cos(x2)*sin(x3)  ,  .  .  . 

(mue /x4~2) ~(l/3)*sin(x2)*cos(x3)  ,  .  .  . 

-2/3/ ( mue / x4  ~ 2) ~(2/3)*sin(x2)*sin(x3) *mue /x4 . .  . 

-3]  ; 

“/Using  Xn  calculate  r.target 
r_targetl  =  r.target (Xn) ; 

Ro  =  [  Rx  ( i  ,  1 )  ;  Ry  ( i  ,  1 )  ;  Rz  ( i  ,  1 )  ]  ; 
rho_eci=r_targetl -Ro ; 

R  =  Ro  ; 

R_transform  =  coordinate.transf orm (R) ; 
rho_sez=R_transform*rho_eci ; 
ps=rho_sez(l  ,1)  ; 
pe  =  rho_sez (2 , 1)  ; 
pz  =  rho_sez (3 , 1 )  ; 

“/Obtain  Az.calc  and  El_calc 
Az=atan2(pe ,-ps) ; 

El=asin(pz/ sqrt (ps~2+pe~2+pz~2) ) ; 

Az_calc (i , 1) =Az ; 

El_calc (i  ,  1) =E1  ; 

“/Calculate  H 

H_tilde=  [pe/ (pe ~2  +  ps ~2)  ,  -ps/ (pe ~2  +  ps ~2)  ,  0;... 

-pz/ (pe~2+ps~2+pz~2) ~ (3/2) *ps  / (l-pz~2/ (pe~2+ps~2+pz.  .  . 
~2) ) “  (1/2)  ,  -pz/(pe"2  +  ps~2  +  pz~2)  “(3/2)*pe/(l-pz... 

“2/ (pe"2  +  ps~2  +  pz~2) ) ~  (1/2)  , (pe~2  +  ps~2)  / (pe~2+ps~2+.  .  . 
pz~2)~(3/2)/((pe"2+ps~2)/(pe~2+ps~2+pz~2))~(l/2)]; 
H=H_tilde*R_transform ; 

“/Calculate  T 
T  =  H*pseudo_phi*phi  ; 
if  i==l 

BigT=T; 

else 

BigT= [BigT ; T] ; 

end  ; 


“/Calculate  residuals 
residual  =  zeros  (2 , 1)  ; 

149  residual (1 , 1) =Az_ob (i , 1) -Az_calc  (i  ,  1)  ; 

residual (2 , 1) =El_ob (i , 1) -El_calc (i  ,  1)  ; 
if  i==l 

BigR=residual ; 

else 

154  BigR= [BigR ; residual] ; 

end 

end 


“/  Update  State  Estimate  (at  epoch.) 
159  P= inv ( BigT ’* Qinv * BigT )  ; 
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delta_X=P*BigT ’*Qinv*BigR; 
Xo=Xo+delta_X . * [0 . 2  0.2  0.2  0.2]’; 


°/0  Check  for  convergence  by  looking  at  percent  change  in  the 
residual  norm 
164  res_norm=norm (BigR) 

res_norm_hist=[res_norm_hist ; res.norm] ; 
if  abs (res_norm-res_norm_old) / res_norm  <0 . 0001 
ANS  =1; 

disp ( 1  converged ’ ) 

169  end 

if  count >max_count 
ANS  =  1  ; 

disp ( 1  Exceeded  Max  Count1) 

174  end 

res_norm_old=res_norm ; 
end 

179  Xo (1 , 1) =Xo (1  ,  1) *180/pi ; 

Xo  (2 , 1)  =Xo  (2 ,1)  *  180/ pi  ; 

Xo (3 , 1) =Xo (3,1) *180/ pi ; 
display ( Xo ) 
display (P) 

184 

'/.The  loops  below  are  seperating  out  azimuth  and  elevation  .. 
residuals 

for  i=l : length (BigR) 
odd= i *2-1 ; 
oddl ( i , 1 ) =  odd ; 

189  even=i*2-2; 

evenl (i ,l)=even; 

end 


for  i=l : length (BigR) 

194  for  a=l : length ( evenl ) 

if  i==evenl (a  ,  1) 
el  =  BigR ( i  ,  1)  ; 
ell (a , 1) =el ; 
end 

199  end 

end 

for  i=l : length (BigR) 

for  a=l : length ( oddl ) 
if  i==oddl (a  ,  1) 

204  Az  =  BigR (i  ,  1)  ; 

Azl (a , 1) =Az ; 
end 
end 


end 


Appendix  C.  Matlab  Functions 

This  appendix  contains  all  function  code  referenced  in  the  core  code  in  appendix  1. 


C.l  Observation  Function 

Listing  C.l:  The  observation  function  G  used  in  NLleastsquares.m. 
(appendix3/G.m) 

V.  ... 


2  / 
7. 
7. 
7. 
7. 

7  1 

% 

% 

7. 

7. 

12  7. 


funct ion  G . m 


this  is  the  observation  funtion  for  least  squares. 


input  s 

rx , ry  ,  rz 
Rx  ,  Ry  ,  Rz 
output  s 
Az_ob 
El_ob 


description 
-  position  of  target 
-position  of  sensor 

Observed  Azimuth 
Observed  Elevation 


range  /  units 
ECI  KM 
ECI  KM 

radians 

radians 


7«  [Az_ob  ,  El_ob  ]  =G  (  day  ,  month  ,  year  ,  hour  ,  sec  ,  rx  ,  ry  ,  rz  ,  Rx  ,  Ry  ,  Rz  ) 
7.  files  needed  to  run:  coordinate.transf on 


17  7o[Az_ob,  El_ob]  =G  (rx  ,  ry  ,  rz  ,  Rx  ,  Ry  ,  Rz) 

7.  ... 


function  [Az_ob  ,  El_ob]  =G  (rx  ,  ry  ,  rz  ,  Rx  ,  Ry  ,  Rz) 

22 


7«Form  vectors  in  ECI 
R=  [Rx  ,  Ry  ,  Rz]  ; 
r= [rx , ry , rz]  ; 

27 

7>rho  in  the  eci 
for  i=l : length (R) 

A  =  r  ( i  ,  : ) -R ( i ,  : )  ; 
rho.eci (i , : ) =A ; 

32  end 

7«rho  in  the  sez  frame 
for  i=l : length (rho.eci ) 

R_ transform  =  coordinate_transform(R(i ,  : )  ’)  ; 
37  B  =  R_transform*rho_eci(i,  :)  ’  ; 

rho.sez (i , : ) =B ; 

end 


87 


'/.generate  az  and  el  data 
42  ps=rho_sez(:,l); 
pe=rho_sez(:  ,2)  ; 
pz=rho_sez(:  ,3)  ; 
for  i  =  1  :  length.  (  rlio_sez  ) 

F=  at an2 (pe(i,l)  ,-ps(i,l)); 

47  Az ( i , 1 ) =F ; 

end 

for  i  =  1  :  length  (  rlio_sez  ) 

B=asin (pz (i , 1) /sqrt (ps (i , 1) ~2  +  pe  (i  ,  1) ~2  +  pz (i , 1) ~2) )  ; 
52  El ( i , 1 ) =B ; 

end 

Az_ob  =  Az  ; 

El _ob  =  E1  ; 

C.2  Target  function 

Listing  C.2:  The  target  function. (appendix3/rtarget.m) 

7.  ... 


2  7. 
7. 
7. 
7. 


7. 

7  7. 
7. 
7. 
7. 
7. 

12  7. 


function  r_target.m 


this  function  finds  the  calculated  position  of  the  target 
Satellite  from  Xn . 


input  s 
Xn 
r 


description 

State  Vector  at  time  n 
height  value  being  used 


range  /  units 

varied 

KM 


output  s 
r _t arget 


target  position  in  ECI  fram  KM 


7.  r.target  =  r.target  (Xn  ,  r  ,  mue) 

7.  ••• 


function  r_targetl  =  r_target(Xn) 

mue=398601 ; 

Omega=Xn (1,1); 

22  I  =  Xn  (2 ,1)  ; 

Theta  =  Xn (3,1)  ; 

Thetadot =Xn  (4 , 1 )  ; 

RIn=  [1  0  0; 

0  cos  (I)  -sin  (I); 

0  sin (I)  cos  (I)]; 
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ROmegan= [cos ( Omega)  -sin(Omega)  0; 
sin(Omega)  cos(Omega)  0; 

0  0  1]; 

R_pqw=[cos (Theta) ; 

32  sin (Theta) ; 

0]  ; 

n= (mue /Thetadot~2) “(1/3)  ; 
r_targetl=ROmegan*RIn*R_pqw*n; 

C.3  Julian  Day 

Listing  C.3:  The  julian  day  function  for  reference. (appendix3/JDAY.m) 

7.  ... 


7. 

7,  function  jday.m 

7. 


5  °/ 

i  this  function 

finds  the  julian 

date  given  the 

year  , 

month  ,  day  , 

°y 

and  time  . 

l 

7 

inputs 

description 

range 

/  units 

7 

1  year 

-  year 

1900 

. .  2100 

7 

1  mon 

-  month 

1  .  . 

12 

10  7 

l  day 

-  day 

1  .  . 

28 ,29,30,31 

7 

l  hr 

-  universal  time 

hour 

0  .  . 

23 

7 

1  min 

-  universal  time 

min 

0  .  . 

59 

7 

1  sec 

-  universal  time 

sec 

0  .  OdO 

..  59.999. 

dO 

/i 

15  7 

i  outputs 

7 

>  j  d 

-  julian  date 

days 

from  4713  . 

be 

7 

l  jd  =  jday (year . 

,  month  ,  day  ,  hour 

,  min ,  sec ) 

7. 


function  JD  =  jday(year,  month,  day,  hour,  min,  sec) 

JD=367*year -floor ( (7* (year  +  floor ( (month+9) /12) ) ) /4) tfloor ( (275* .  .  . 
mo nth)/9)+day +1721013. 5+((((sec/60)+min)/60)+hour)/24; 

C-4  Greenwich  Mean  Sidereal  Time 

Listing  C.4:  The  Ogmst  function  ,(appendix3/GSTIME.m) 

7.  ... 


7. 


function  gstime 


7. 

4  °/0  this  function  finds  the  greenwich  sidereal  time. 

7. 


% 

input  s 

description 

range 

/ 

units 

7. 

j  dut  1 

-  julian  date  of  utl 

days 

from  4713 

7 

be 

/o 

9  7. 

output  s 

7. 

y 

gst 

-  greenwich  sidereal  time 

0  to 

2pi 

rad 

/o 

7. 

locals 

7. 

temp 

-  temporary  variable  for 

reals 

rad 

14  7. 

tut  1 

-  julian  centuries  from  the 

7. 

j  an  1,  2000  12  h  epoch 

(utl  ) 

% 

ref  erences 

7. 

vallado 

2001,  191,  eq  3-45 

7,  gst  =  gstime(JD); 
19  7.  ... 


function  gst  =  gstime(JD); 

twopi  =  2 . 0*pi ; 

deg2rad  =  pi/180.0; 

7,  -  implementation 


tutl=  (  JD  -  2451545.0  )  /  36525.0; 


temp  =  -  6.2e-6  *  tutl  *  tutl  *  tutl  +  0.093104  *  tutl  * 
tut  1  ... 

29  +  (876600.0  *  3600.0  +  8640184.812866)  *  tutl  +  . 

67310 . 54841 ; 

temp  =  modulo (  t emp * deg2r ad /240 . 0 , twopi  ); 

7,  -  check  quadrants  ... 


34  if  (  temp  <  0.0  ) 

temp  =  temp  +  twopi ; 
end 


gst  =  temp ; 


C.  5  Coordinate  Transform 

Listing  C.5:  The  coordinate  transform  function. 
(appendix3  /  coordinatetransform.m) 

17.  ••• 


7. 
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7. 

% 

7. 

6  % 
7. 
7. 
7. 
% 

11  7. 


7. 

7. 


function  coordinate.transf orm.m 

this  function  finds  the  coordinate  transform  from  ECI  to  SEZ . 

inputs  description  range  /  units 

R  position  vector  of  sensor  KM  ECI 

outputs  : 

R.transform  coordinate  transfrom  unitless 

R.transf orm  =  coordinate.transf orm (R) 


16  function  R.transform  =  coordinate.transf orm (R) 

‘/.Coordinate  Transform 
z=R/ norm ( R) ; 
kvec  =  [0 ; 0  ;  1]  ; 

21  b=cross (kvec ,R) ; 
bl=norm (b) ; 
e=b/bl ; 
s  =  cross (e  ,  z)  ; 

C=  [s  ’  ;  e  ’  ;  z  ’  ]  ; 

26  R_transf orm=C ; 

C.  6  Modulus 

Listing  C.6:  Simple  but  necessary  modulus  function.  (appendix3/modulo.m) 

%  ... 


2  7. 

7«  modulo  f  unct  ion 

7. 


7. 


y. 

7  % 

y. 

7. 

7. 

7. 

12  7. 

% 


this  function  finds  the  modulous  value.  The 
or  negative . 


input  s 
x 

xv 


description 

-  input  argument 

-  input  argument  to  mod  with 


output  s 

y 

y  =  modulo 


-  answer 
(x , xv )  ; 


7. 


result  is  positive 


range  /  units 


function  y  =  modulo  (x,xv); 

17  y  =  x  -  xv  *  fix(x  /  xv)  ; 
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C.  7  Search  Algorithm 


Listing  C.7:  The  search  algorithm  uses  a  modified  single  line  version  of  the  catalog. 
(appendix3/simplesearch.m) 
clear  all 
2  clc 

format  long 

*/,f  id  =  f  open  (  1  2ndlinecatalog  .  txt  ’  ,  ’r’)  ;'/,the  one  without  debris 

fid  =  f  open  (  ’  catalog2ndlinewdebr  is  .  txt  1  ,  ’r ’)  ;'/,the  one  with  debris 
7  dataimport  =  fscanf(fid,  ’7,  g  7.g  7«g  7,g  7.g  7.g  7.g  7.g  ’  ,  [8  inf]); 

dataimport=dataimport ’ ; %  It  has  two  rows  now. 
fclose (fid) 
line=dataimport (: ,1) ; 
de si gnat or  =  dataimport ( :  ,2)  ; 

12  inclination=dataimport ( : ,3) ; 
raan  =  dataimport ( :  ,4)  ; 
eccentric ity  =  dat a import ( :  ,5)  ; 
argument _of_perigee=dataimport ( : ,6) ; 
mean.anomaly  =  dataimport  (  :  ,7)  ; 

17  mean.mot ion  =  dataimport ( :  ,8)  ; 

fid  =  f open (’ 2ndlinecatalog_epoch . txt  1  ,  ’r’); 

dataimport  =  fscanf(fid,  ’  7,14g>,  [1  inf]); 

dat  aimport  =  dat  aimport  ’  ;  %  It  has  two  rows  now. 

22  fclose (fid) 

Epoch  =  d at aimport (:  ,1)  ; 
f igur e 
hold  on 

plot (inclination , raan ,  ’  .  ’  ) 

27  title (’ Inclination  vs  RAAN’) 
xlabel ( ’ Inclination ’ ) 
ylabel ( ’ RAAN ’ ) 
set (gca ,  ’ XTick  ’  ,0:10:150) 
hold  off 
32 

f  ound  =0 ; 

for  i = 1 : length ( 1 ine ) 

if  inclination (i  ,  1)  >  98.5  kk  inclination  (i  ,  1)  <  99.5  kk  raan(... 

i,l)  >  65.4  kk  raan(i,l)  <  66.4 
37  f ound  =  f ound  + 1  ; 

match(found)=i ; 

end 

end 

if  found==0,  disp(’No  mat che s 1 ) , end 
42  if  found~=0 

disp([;Found  num2str (found) , ’  matches.’]) 
disp (  [  1  ’  ]  ) 

for  i=l : found 

disp ([’ Satellite  Number:  num2str (designator (match (1 , i ..  . 

)  ,1))]) 
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47 


,  num2str ( inclination (match (1 , 


52 


disp ( [ ’ Inclination : 
i)  ,1))]) 

disp([’RAAN:  ’  ,  num2str (raan(match(l,i) ,1)). 

]) 

disp ([ 1  Eccentricity :  ’,  num2str ( eccentricity (match .  . 

(  1  ,  i )  ,1))]) 

di sp ([’ Argument  of  Perigee:’,  num2str ( argument _ of _perigee ( 
match (1 , i)  ,  1) ) ] ) 

disp ([’Mean  Anomaly:  ’,  num2str (mean.anomaly (match .  . 

(  1  ,  i )  ,1))]) 

disp([’Mean  Motion:  ’,  num2str (mean.motion (match (1 , 

i)  ,1)  ,  ’*/.2.8f  ’)]) 
disp  (  [  ’  ’  ]  ) 

end 

end 


93 


Bibliography 


1.  Briggs,  RE  and  JW  Slowley.  “An  Iterative  Method  of  Orbit  Determination  from 
Three  Observations  of  a  Nearby  Satellite”.  SAO  Special  Report ,  27(part  1),  1959. 

2.  Burns,  R.E.  “Solution  of  the  Angles-Only  Satellite  Tracking  Problem”.  NASA 
Technical  Paper ,  3667:25,  1997. 

3.  Escobal,  P.R.  Methods  of  Orbit  Determination.  Wiley  New  York,  1965. 

4.  Gooding,  RH.  “A  New  Procedure  for  Orbit  Determination  Based  on  Three  Lines 
of  Sight  (Angles  Only)”.  DRA/RAE  Technical  Report ,  93004,  1993. 

5.  Gooding,  RH.  “A  new  procedure  for  the  solution  of  the  classical  problem  of 
minimal  orbit  determination  from  three  lines  of  sight”.  Celestial  Mechanics  and 
Dynamical  Astronomy ,  66(4) :387— 423,  1996. 

6.  Kristensen,  LK.  “Orbit  Determination  by  Four  Observations”.  ASTRONOMIS- 
CHE  NACHRICHTEN ,  316(4):261,  1995. 

7.  Li,  Q.,  F.  Guo,  J.  Li,  and  Y.  Zhou.  “Research  of  Satellite-to-Satellite  Passive 
Tracking  Using  Bearings-Only  Measurements  In  J2000  ECI  Frame”.  Radar,  2006. 
CIE’06.  International  Conference  on,  1-4,  2006. 

8.  Marsden,  BG.  “Initial  orbit  determination- The  pragmatist’s  point  of  view”.  The 
Astronomical  Journal,  90:1541,  1985. 

9.  OSEDACZ,  R.  “Orbit  determination  of  sunlight  illuminated  objects  detected  by 
overhead  platforms(M.  S.  Thesis)”.  1989. 

10.  Prussing,  J.E.  and  B.A.  Conway.  Orbital  Mechanics.  Oxford  University  Press 
New  York,  1993. 

11.  SABOL,  C.,  K.T.  ALFRIEND,  and  D.  WIESE.  “Angles-only  orbit  updates  for 
low-earth-orbit  satellites” .  Journal  of  guidance,  control,  and  dynamics ,  30(2):314- 
321,  2007. 

12.  Sabol,  C.  and  R.  Culp.  “Improved  angular  observations  in  geosynchronous  orbit 
determination”.  Journal  of  Guidance,  Control,  and  Dynamics,  24(1):123-130, 
2001. 

13.  Taff,  LG.  “On  initial  orbit  determination” .  Astronomical  Journal,  89:1426-1428, 
1984. 

14.  Tapley,  B.,  B.  Schutz,  and  G.  Born.  Statistical  Orbit  Determination.  Elsevier 
Academic  Press,  Boston,  2004. 

15.  Vallado,  D.A.  Fundamentals  of  Astrodynamics  and  Applications.  Kluwer  Aca¬ 
demic  Pub,  2001. 


94 


16.  Wiesel,  William  E.  Spaceflight  Dynamics.  McGraw-Hill,  Boston,  1989. 

17.  Wiesel,  William  E.  Modern  Methods  of  Orbit  Determination.  Apahelion  Press, 
BeaverCreek  Ohio,  2003. 


95 


Index 


The  index  is  conceptual  and  does  not  designate  every  occurrence  of  a  key¬ 
word.  Page  numbers  in  bold  represent  concept  definition  or  introduction. 


Code  Listings 

appendix2/GaussIOD.m,  71 
appendix2/NLleastsquares.m,  82 
appendix2/SphericaltrigIODm.m,  78 
appendix3/G.m,  87 
appendix3/GSTIME.m,  89 
appendix3/JDAY.m,  89 
appendix3/coordinatetransform.m,  90 
appendix3/modulo.m,  91 
appendix3/rtarget.m,  88 
appendix3/simplesearch.m,  92 


96 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704—0188),  1215  Jefferson  Davis  Highway, 
Suite  1204,  Arlington,  VA  22202—4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection 
of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 

10-03-2008 


2.  REPORT  TYPE 

Master’s  Thesis 


3.  DATES  COVERED  (From  —  To) 

Sept  2007  Mar  2008 


4.  TITLE  AND  SUBTITLE 


Space  Based  Satellite  Tracking  and  Characterization 
Utilizing  Non-Imaging  Passive  Sensors 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


Bradley  R.  Townsend,  CPT,  USA 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Institute  of  Technology 
Graduate  School  of  Engineering  and  Management 
2950  Hobson  Way 
WPAFB  OH  45433-7765 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


AFIT/G  A /ENY/08-M06 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approval  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

A  technique  is  developed  to  determine  the  orbit  of  a  sunlight  illuminated  satellite  passing  through  the  field-of-view  of  a 
sensor  platform  in  a  Highly  Elliptical  Orbit  (HEO)  and  Geosynchronous  orbit  (GEO).  The  technique  develops  two 
different  methods  of  initial  orbit  determination.  The  first  relies  on  the  Gauss  initial  orbit  determination  method  to 
develop  an  estimate  of  the  state  from  angular  data.  The  second  method  relies  on  positional  data  of  the  target  relative  to 
the  Earth’s  background  to  determine  an  estimate  of  the  state.  These  estimates  are  then  refined  in  a  non-linear  least 
squares  routine.  This  estimate  of  the  state  is  then  used  to  identify  the  target  from  the  Air  Force  Space  Command 
satellite  catalog. 


15.  SUBJECT  TERMS 

Gauss,  fastwalker,  Angles  only,  Orbit  Determination,  Highly  Elliptical  Orbit,  Geosynchronous  orbit 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF 

PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

LTC  Nathan  Titus 

u 

u 

u 

uu 

111 

19b.  TELEPHONE  NUMBER  (include  area  code) 

(937)  785-3636,  ext  7469 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


